From bb607377b1b4d4bbdb44f8ee61e34c35f3a8391f Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Thu, 11 Apr 2024 18:49:18 -0400 Subject: [PATCH 01/13] Optimize and slightly refactor sorteachterm. --- src/treetensornetworks/opsum_to_ttn.jl | 38 +++++++++++++++++--------- 1 file changed, 25 insertions(+), 13 deletions(-) diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl index b2c2a1f0..0ac80858 100644 --- a/src/treetensornetworks/opsum_to_ttn.jl +++ b/src/treetensornetworks/opsum_to_ttn.jl @@ -6,14 +6,7 @@ using ITensors.ITensorMPS: ITensorMPS # 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} @@ -64,9 +57,9 @@ function ttn_svd( 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 + vs = _default_vertex_ordering(sites,root_vertex) # ToDo: Add check in ttn_svd that the ordering matches that of find_index_in_tree, which is used in sorteachterm #fermion-sign! - es = reverse(reverse.(post_order_dfs_edges(sites, root_vertex))) # store edges in fixed ordering relative to root + es = _default_edge_ordering(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 @@ -402,11 +395,29 @@ function computeSiteProd(sites::IndsNetwork{V,<:Index}, ops::Prod{Op})::ITensor return T end +function _default_vertex_ordering(g::AbstractGraph,root_vertex) + return reverse(post_order_dfs_vertices(g, root_vertex)) +end + +function _default_edge_ordering(g::AbstractGraph,root_vertex) + return reverse(reverse.(post_order_dfs_edges(g, root_vertex))) +end + # This is almost an exact copy of ITensors/src/opsum_to_mpo_generic:sorteachterm except for the site ordering being # given via find_index_in_tree # changed `isless_site` to use tree vertex ordering relative to root function sorteachterm(os::OpSum, sites::IndsNetwork{V,<:Index}, root_vertex::V) where {V} + println("in sort") os = copy(os) + ordering = _default_vertex_ordering(sites,root_vertex) + site2pos=Dict(zip(ordering,1:length(ordering))) + # 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) + return site2pos[site] + end + function find_index_in_tree(o::Op, g::AbstractGraph, root_vertex) + return find_index_in_tree(ITensors.site(o), g, root_vertex) + end findpos(op::Op) = find_index_in_tree(op, sites, root_vertex) isless_site(o1::Op, o2::Op) = findpos(o1) < findpos(o2) N = nv(sites) @@ -468,6 +479,7 @@ function sorteachterm(os::OpSum, sites::IndsNetwork{V,<:Index}, root_vertex::V) t *= ITensors.parity_sign(perm) ITensors.terms(os)[n] = t end + println("out sort") return os end @@ -490,10 +502,10 @@ function TTN( is_leaf(sites, root_vertex) || error("Tree root must be a leaf vertex.") os = deepcopy(os) - os = sorteachterm(os, sites, root_vertex) - os = ITensorMPS.sortmergeterms(os) # not exported + @time os = sorteachterm(os, sites, root_vertex) + @time os = ITensorMPS.sortmergeterms(os) # not exported if algorithm == "svd" - T = ttn_svd(os, sites, root_vertex; kwargs...) + @time T = ttn_svd(os, sites, root_vertex; kwargs...) else error("Currently only SVD is supported as TTN constructor backend.") end From 6884310c1ab35d01151839cecd3e0df7c06f0b2e Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 12 Apr 2024 13:13:17 -0400 Subject: [PATCH 02/13] Optimized the graph logic in ttn_svd further. --- src/treetensornetworks/opsum_to_ttn.jl | 329 +++++++++++++++++-------- 1 file changed, 225 insertions(+), 104 deletions(-) diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl index c2ee2756..48926722 100644 --- a/src/treetensornetworks/opsum_to_ttn.jl +++ b/src/treetensornetworks/opsum_to_ttn.jl @@ -34,12 +34,25 @@ using StaticArrays: MVector # determine 'support' of product operator on tree graph function span(t::Scaled{C,Prod{Op}}, g::AbstractGraph) where {C} - spn = eltype(g)[] + spn = Set{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])) + nterms == 1 && return [ITensors.site(t[1])] + #sites=[ITensors.site(t[i]) for i in 1:nterms] + # ToDo: figure out better implementation, maybe via steiner_tree + + #steiner_spn=ITensorNetworks.NamedGraphs.steiner_tree(g,sites) + for i in 1:nterms, j in i+1:nterms + path = Set(vertex_path(g, ITensors.site(t[i]), ITensors.site(t[j]))) spn = union(spn, path) end + #nspn = eltype(g)[] + #for i in 1:nterms, j in i:nterms + # path = vertex_path(g, ITensors.site(t[i]), ITensors.site(t[j])) + # nspn = union(spn, path) + #end + #@assert nspn==spn + #@show vertices(steiner_spn), spn + #@assert issetequal(vertices(steiner_spn),spn) return spn end @@ -115,6 +128,8 @@ function ttn_svd( for v in vs is_internal[v] = isempty(sites[v]) if isempty(sites[v]) + # FIXME: This logic only works for trivial flux, breaks for nonzero flux + # ToDo: add assert or fix and add test! sites[v] = [Index(Hflux => 1)] end end @@ -124,113 +139,216 @@ function ttn_svd( 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) - + #ToDo: precompute span of each term and store + # compute span of each term + spans=Dict{eltype(os),Set{vertextype_sites}}() + @time "span" begin + for term in os + spans[term]=span(term,sites) + end + end # 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 = calc_qn(incoming) - not_incoming_qn = calc_qn(not_incoming) - 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{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 - if !isempty(incoming) - # get the correct map from edge=>QN to term and channel - # this checks if term exists on edge=>QN ( otherwise insert it) and returns it's index - 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 = ITensorMPS.posInLink!(cinmap, incoming) - bond_col = ITensorMPS.posInLink!(coutmap, not_incoming) # get incoming channel - bond_coef = convert(coefficient_type, ITensors.coefficient(term)) - q_inbond_coefs = get!( - inbond_coefs[edge_in], incoming_qn, ITensorMPS.MatElem{coefficient_type}[] - ) - push!(q_inbond_coefs, ITensorMPS.MatElem(bond_row, bond_col, bond_coef)) - T_inds[dim_in] = bond_col - T_qns[dim_in] = -incoming_qn + @time "first" begin + for v in vs + #copy(sites) + _g=deepcopy(ITensorNetworks.underlying_graph(sites)) + NamedGraphs.rem_vertex!(_g,v) + subgraphs=NamedGraphs.connected_components(_g) + # for every vertex, find all edges that contain this vertex + # ToDo: use neighborhood instead of going through all edges + # slight complication is to preserve directionality of edge? + alledges=es + edges=edgetype_sites[] + edge_ordering=Int[] + + for other_v in NamedGraphs.neighbors(sites,v) + e=edgetype_sites(v,other_v) + if e in es + push!(edges,e) + push!(edge_ordering, findfirst(isequal(e),es)) + else + @assert reverse(e) in es + push!(edges,reverse(e)) + push!(edge_ordering,findfirst(isequal(reverse(e)),es)) + end end - for dout in dims_out - coutmap = get!( - outmaps, edges[dout] => outgoing_qns[edges[dout]], Dict{Vector{Op},Int}() - ) - T_inds[dout] = ITensorMPS.posInLink!(coutmap, outgoing[edges[dout]]) # add outgoing channel - T_qns[dout] = outgoing_qns[edges[dout]] + edges=edges[sortperm(edge_ordering)] + #edges=sort(edges,by=index(es)) + #ref_edges = filter(e -> dst(e) == v || src(e) == v, es) + #@assert issetequal(edges,ref_edges) + # 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] + if !(typeof(edge_in)<:Vector) + v_in=only(setdiff([src(edge_in),dst(edge_in)],[v,])) end - # if term starts at this site, add its coefficient as a site factor - site_coef = one(coefficient_type) - if (isnothing(dim_in) || T_inds[dim_in] == -1) && - ITensors.argument(term) ∉ site_coef_done - site_coef = ITensors.coefficient(term) - site_coef = convert(coefficient_type, site_coef) # required since ITensors.coefficient seems to return ComplexF64 even if coefficient_type is determined to be real - push!(site_coef_done, ITensors.argument(term)) + v_outs=vertextype_sites[] + for e_out in edges_out + push!(v_outs,only(setdiff([src(e_out),dst(e_out)],[v,]))) 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!") + in_graphs=[] + out_graphs=[] + which_edges=edgetype_sites[] + _vs=vertextype_sites[] + for subgraph in subgraphs + is_incoming=false + typeof(edge_in)<:Vector || (is_incoming = v_in in subgraph) + if !is_incoming + + for i in eachindex(v_outs) + if v_outs[i] in subgraph + _e = edges_out[i] + break + end + end else - push!(onsite, Op("Id", v)) + _e=edge_in end + append!(_vs,subgraph) + append!(which_edges,fill(_e,length(subgraph))) end - # save indices and value of symbolic tensor entry - el = QNArrElem(T_qns, T_inds, site_coef * Prod(onsite)) + + which_subtree=Dict{vertextype_sites,edgetype_sites}(zip(_vs,which_edges)) + + # 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 + v in spans[term] || continue + #crosses_vertex(term, sites, v) || continue + + # filter out factors that come in from the direction of the incoming edge + # ToDo: edge_path could compute over spanning tree of term instead of full graph + # not sure if this affects the performance of edge_path + # ToDo: possibly precompute edge_path once for all entrie in term + # and reuse in each of the calls + #@show term + #@show v + #edge_paths=map(t->edge_path(sites, ITensors.site(t), v),ITensors.terms(term)) + #in principle converting to set should be faster since O(1) access + #edge_paths=map(t->Set(edge_path(sites, ITensors.site(t), v)),ITensors.terms(term)) + + #@show typeof(edge_paths) + #for i in 1:length(ITensors.terms(term)) + # @show i,ITensors.terms(term)[i] + #end + factors=ITensors.terms(term) + onsite = filter(t -> (ITensors.site(t) == v), factors) + not_onsite_factors=setdiff(factors,onsite) + #incoming= ITensors.terms(term)[map(_edges -> edge_in in _edges,edge_paths)] + incoming=filter(t->which_subtree[ITensors.site(t)]==edge_in, not_onsite_factors) + #@show incoming + #@show incoming_alt + #@assert incoming==incoming_alt + #incoming = filter( + # t -> edge_in ∈ edge_path(sites, ITensors.site(t), v), ITensors.terms(term) + #) + + #incoming_alt = ITensors.terms(term)[(in.(edge_in,edge_paths))...] + #@show incoming + #@show incoming_alt + #@assert incoming_alt==collect(incoming) + # also store all non-incoming factors in standard order, used for channel merging + #show v + #not_incoming = ITensors.terms(term)[map(_edges -> edge_in ∉ _edges,edge_paths)] + + not_incoming = filter(t-> (ITensors.site(t) == v) || which_subtree[ITensors.site(t)]!=edge_in,factors) + #not_incoming = filter( + # t -> edge_in ∉ edge_path(sites, ITensors.site(t), v), ITensors.terms(term) + #) + # filter out factor that acts on current vertex + + + # 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 + #) + #outgoing = Dict( + # e => ITensors.terms(term)[map(_edges -> e in reverse.(_edges), edge_paths)] + # for e in edges_out + #) + outgoing = Dict( + e => filter(t->which_subtree[ITensors.site(t)]==e,not_onsite_factors) + for e in edges_out + ) + #@assert outgoing==outgoing_alt + # compute QNs + incoming_qn = calc_qn(incoming) + not_incoming_qn = calc_qn(not_incoming) + 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{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 + if !isempty(incoming) + # get the correct map from edge=>QN to term and channel + # this checks if term exists on edge=>QN ( otherwise insert it) and returns it's index + 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 = ITensorMPS.posInLink!(cinmap, incoming) + bond_col = ITensorMPS.posInLink!(coutmap, not_incoming) # get incoming channel + bond_coef = convert(coefficient_type, ITensors.coefficient(term)) + q_inbond_coefs = get!( + inbond_coefs[edge_in], incoming_qn, ITensorMPS.MatElem{coefficient_type}[] + ) + push!(q_inbond_coefs, ITensorMPS.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 + coutmap = get!( + outmaps, edges[dout] => outgoing_qns[edges[dout]], Dict{Vector{Op},Int}() + ) + T_inds[dout] = ITensorMPS.posInLink!(coutmap, outgoing[edges[dout]]) # add outgoing channel + T_qns[dout] = outgoing_qns[edges[dout]] + end + # if term starts at this site, add its coefficient as a site factor + site_coef = one(coefficient_type) + if (isnothing(dim_in) || T_inds[dim_in] == -1) && + ITensors.argument(term) ∉ site_coef_done + site_coef = ITensors.coefficient(term) + site_coef = convert(coefficient_type, site_coef) # required since ITensors.coefficient seems to return ComplexF64 even if coefficient_type is determined to be real + 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 + push!(tempTTN[v], el) + end - ITensorMPS.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 = ITensorMPS.toMatrix(mat) - U, S, V = svd(M) - P = S .^ 2 - truncate!(P; maxdim, cutoff, mindim) - tdim = length(P) - nc = size(M, 2) - Vs[edges[dim_in]][q] = Matrix{coefficient_type}(V[1:nc, 1:tdim]) + ITensorMPS.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 = ITensorMPS.toMatrix(mat) + U, S, V = svd(M) + P = S .^ 2 + truncate!(P; maxdim, cutoff, mindim) + tdim = length(P) + nc = size(M, 2) + Vs[edges[dim_in]][q] = Matrix{coefficient_type}(V[1:nc, 1:tdim]) + end end end end - # compress this tempTTN representation into dense form link_space = Dict{edgetype_sites,Index}() @@ -256,6 +374,7 @@ function ttn_svd( for v in vs # redo the whole thing like before + # ToDo: use neighborhood instead of going through all edges, see above 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) @@ -350,9 +469,10 @@ function ttn_svd( if is_internal[v] H[v] += iT else - if hasqns(iT) - @assert flux(iT * Op) == Hflux - end + #ToDo: Remove this assert since it seems to be costly + #if hasqns(iT) + # @assert flux(iT * Op) == Hflux + #end H[v] += (iT * Op) end end @@ -433,11 +553,12 @@ end function sorteachterm(os::OpSum, sites::IndsNetwork{V,<:Index}, root_vertex::V) where {V} println("in sort") os = copy(os) - ordering = _default_vertex_ordering(sites,root_vertex) - site2pos=Dict(zip(ordering,1:length(ordering))) + # linear ordering of vertices in tree graph relative to chosen root, chosen outward from root + ordering = _default_vertex_ordering(sites,root_vertex) + site_positions=Dict(zip(ordering,1:length(ordering))) function find_index_in_tree(site, g::AbstractGraph, root_vertex) - return site2pos[site] + return site_positions[site] end function find_index_in_tree(o::Op, g::AbstractGraph, root_vertex) return find_index_in_tree(ITensors.site(o), g, root_vertex) From c7fa14158cb1b973df4bc34a80c5c194c1a6ab66 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 12 Apr 2024 14:43:33 -0400 Subject: [PATCH 03/13] Cleanup. --- src/treetensornetworks/opsum_to_ttn.jl | 180 +++++++------------------ 1 file changed, 52 insertions(+), 128 deletions(-) diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl index 7c8d48e0..b587caa1 100644 --- a/src/treetensornetworks/opsum_to_ttn.jl +++ b/src/treetensornetworks/opsum_to_ttn.jl @@ -17,25 +17,13 @@ using StaticArrays: MVector # determine 'support' of product operator on tree graph function span(t::Scaled{C,Prod{Op}}, g::AbstractGraph) where {C} - spn = Set{eltype(g)}() + spn = eltype(g)[] nterms = length(t) nterms == 1 && return [ITensors.site(t[1])] - #sites=[ITensors.site(t[i]) for i in 1:nterms] - # ToDo: figure out better implementation, maybe via steiner_tree - - #steiner_spn=ITensorNetworks.NamedGraphs.steiner_tree(g,sites) for i in 1:nterms, j in i+1:nterms - path = Set(vertex_path(g, ITensors.site(t[i]), ITensors.site(t[j]))) + path = vertex_path(g, ITensors.site(t[i]), ITensors.site(t[j])) spn = union(spn, path) end - #nspn = eltype(g)[] - #for i in 1:nterms, j in i:nterms - # path = vertex_path(g, ITensors.site(t[i]), ITensors.site(t[j])) - # nspn = union(spn, path) - #end - #@assert nspn==spn - #@show vertices(steiner_spn), spn - #@assert issetequal(vertices(steiner_spn),spn) return spn end @@ -44,6 +32,40 @@ function crosses_vertex(t::Scaled{C,Prod{Op}}, g::AbstractGraph, v) where {C} return v ∈ span(t, g) end +# map all vertices `w` of `g` except for `v` to the incident edge of `v` +# which lies in edge_path(g,w,v) +function map_vertices_to_incident_edges(g::AbstractGraph,v,es) + _g=deepcopy(ITensorNetworks.underlying_graph(g)) + NamedGraphs.rem_vertex!(_g,v) + subgraphs=NamedGraphs.connected_components(_g) + + vs=vertextype(g)[] + for e in es + push!(vs,only(setdiff([src(e),dst(e)],[v,]))) + end + _vs=vertextype(g)[] + _es=edgetype(g)[] + + for (e,v) in zip(es,vs) + for i in eachindex(subgraphs) + if v in subgraphs[i] + append!(_vs,subgraphs[i]) + append!(_es,fill(e,length(subgraphs[i]))) + deleteat!(subgraphs,i) + break + end + end + end + @assert isempty(subgraphs) + return Dict(zip(_vs,_es)) +end + +function ordered_incident_edges(g::AbstractGraph,v,ordered_edges) + edges=ITensorNetworks.incident_edges(g,v) + edges=map(e -> e in ordered_edges ? e : reverse(e), edges) + edge_positions=map(e -> findfirst(isequal(e),ordered_edges),edges) + return edges[sortperm(edge_positions)] +end # # Tree adaptations of functionalities in ITensors.jl/src/physics/autompo/opsum_to_mpo.jl # @@ -124,7 +146,7 @@ function ttn_svd( tempTTN = Dict(v => QNArrElem{Scaled{coefficient_type,Prod{Op}},degrees[v]}[] for v in vs) #ToDo: precompute span of each term and store # compute span of each term - spans=Dict{eltype(os),Set{vertextype_sites}}() + spans=Dict{eltype(os),Vector{vertextype_sites}}() @time "span" begin for term in os spans[term]=span(term,sites) @@ -133,133 +155,43 @@ function ttn_svd( # build compressed finite state machine representation @time "first" begin for v in vs - #copy(sites) - _g=deepcopy(ITensorNetworks.underlying_graph(sites)) - NamedGraphs.rem_vertex!(_g,v) - subgraphs=NamedGraphs.connected_components(_g) # for every vertex, find all edges that contain this vertex - # ToDo: use neighborhood instead of going through all edges - # slight complication is to preserve directionality of edge? - alledges=es - edges=edgetype_sites[] - edge_ordering=Int[] - - for other_v in NamedGraphs.neighbors(sites,v) - e=edgetype_sites(v,other_v) - if e in es - push!(edges,e) - push!(edge_ordering, findfirst(isequal(e),es)) - else - @assert reverse(e) in es - push!(edges,reverse(e)) - push!(edge_ordering,findfirst(isequal(reverse(e)),es)) - end - end - edges=edges[sortperm(edge_ordering)] - #edges=sort(edges,by=index(es)) - #ref_edges = filter(e -> dst(e) == v || src(e) == v, es) - #@assert issetequal(edges,ref_edges) + edges=ordered_incident_edges(sites,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] - if !(typeof(edge_in)<:Vector) - v_in=only(setdiff([src(edge_in),dst(edge_in)],[v,])) - end - v_outs=vertextype_sites[] - for e_out in edges_out - push!(v_outs,only(setdiff([src(e_out),dst(e_out)],[v,]))) - end - in_graphs=[] - out_graphs=[] - which_edges=edgetype_sites[] - _vs=vertextype_sites[] - for subgraph in subgraphs - is_incoming=false - typeof(edge_in)<:Vector || (is_incoming = v_in in subgraph) - if !is_incoming - - for i in eachindex(v_outs) - if v_outs[i] in subgraph - _e = edges_out[i] - break - end - end - else - _e=edge_in - end - append!(_vs,subgraph) - append!(which_edges,fill(_e,length(subgraph))) - end - - which_subtree=Dict{vertextype_sites,edgetype_sites}(zip(_vs,which_edges)) + which_incident_edge=map_vertices_to_incident_edges(sites,v,edges) + factor_to_edge(t)=which_incident_edge[ITensors.site(t)] # 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 - v in spans[term] || continue - #crosses_vertex(term, sites, v) || continue - - # filter out factors that come in from the direction of the incoming edge - # ToDo: edge_path could compute over spanning tree of term instead of full graph - # not sure if this affects the performance of edge_path - # ToDo: possibly precompute edge_path once for all entrie in term - # and reuse in each of the calls - #@show term - #@show v - #edge_paths=map(t->edge_path(sites, ITensors.site(t), v),ITensors.terms(term)) - #in principle converting to set should be faster since O(1) access - #edge_paths=map(t->Set(edge_path(sites, ITensors.site(t), v)),ITensors.terms(term)) - #@show typeof(edge_paths) - #for i in 1:length(ITensors.terms(term)) - # @show i,ITensors.terms(term)[i] - #end + v in spans[term] || continue factors=ITensors.terms(term) + + # filter out factor that acts on current vertex onsite = filter(t -> (ITensors.site(t) == v), factors) not_onsite_factors=setdiff(factors,onsite) - #incoming= ITensors.terms(term)[map(_edges -> edge_in in _edges,edge_paths)] - incoming=filter(t->which_subtree[ITensors.site(t)]==edge_in, not_onsite_factors) - #@show incoming - #@show incoming_alt - #@assert incoming==incoming_alt - #incoming = filter( - # t -> edge_in ∈ edge_path(sites, ITensors.site(t), v), ITensors.terms(term) - #) + + # filter out factors that come in from the direction of the incoming edge + incoming=filter(t->factor_to_edge(t)==edge_in, not_onsite_factors) - #incoming_alt = ITensors.terms(term)[(in.(edge_in,edge_paths))...] - #@show incoming - #@show incoming_alt - #@assert incoming_alt==collect(incoming) # also store all non-incoming factors in standard order, used for channel merging - #show v - #not_incoming = ITensors.terms(term)[map(_edges -> edge_in ∉ _edges,edge_paths)] - - not_incoming = filter(t-> (ITensors.site(t) == v) || which_subtree[ITensors.site(t)]!=edge_in,factors) - #not_incoming = filter( - # t -> edge_in ∉ edge_path(sites, ITensors.site(t), v), ITensors.terms(term) - #) - # filter out factor that acts on current vertex - + not_incoming = filter(t-> (ITensors.site(t) == v) || factor_to_edge(t)!=edge_in,factors) # 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 - #) - #outgoing = Dict( - # e => ITensors.terms(term)[map(_edges -> e in reverse.(_edges), edge_paths)] - # for e in edges_out - #) outgoing = Dict( - e => filter(t->which_subtree[ITensors.site(t)]==e,not_onsite_factors) + e => filter(t->factor_to_edge(t)==e,not_onsite_factors) for e in edges_out ) - #@assert outgoing==outgoing_alt + # compute QNs incoming_qn = calc_qn(incoming) not_incoming_qn = calc_qn(not_incoming) @@ -358,7 +290,7 @@ function ttn_svd( for v in vs # redo the whole thing like before # ToDo: use neighborhood instead of going through all edges, see above - edges = filter(e -> dst(e) == v || src(e) == v, es) + edges=ordered_incident_edges(sites,v,es) dim_in = findfirst(e -> dst(e) == v, edges) dims_out = findall(e -> src(e) == v, edges) # slice isometries at this vertex @@ -534,19 +466,12 @@ end # given via find_index_in_tree # changed `isless_site` to use tree vertex ordering relative to root function sorteachterm(os::OpSum, sites::IndsNetwork{V,<:Index}, root_vertex::V) where {V} - println("in sort") os = copy(os) # linear ordering of vertices in tree graph relative to chosen root, chosen outward from root ordering = _default_vertex_ordering(sites,root_vertex) site_positions=Dict(zip(ordering,1:length(ordering))) - function find_index_in_tree(site, g::AbstractGraph, root_vertex) - return site_positions[site] - end - function find_index_in_tree(o::Op, g::AbstractGraph, root_vertex) - return find_index_in_tree(ITensors.site(o), g, root_vertex) - end - findpos(op::Op) = find_index_in_tree(op, sites, root_vertex) + findpos(op::Op) = site_positions[ITensors.site(op)] isless_site(o1::Op, o2::Op) = findpos(o1) < findpos(o2) N = nv(sites) for n in eachindex(os) @@ -607,7 +532,6 @@ function sorteachterm(os::OpSum, sites::IndsNetwork{V,<:Index}, root_vertex::V) t *= ITensors.parity_sign(perm) ITensors.terms(os)[n] = t end - println("out sort") return os end From f7c8507fd04e93ab11e0a101bbd9d5a149f6323c Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 12 Apr 2024 14:50:39 -0400 Subject: [PATCH 04/13] Remove instrumentation. --- src/treetensornetworks/opsum_to_ttn.jl | 210 ++++++++++++------------- 1 file changed, 104 insertions(+), 106 deletions(-) diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl index b587caa1..7e893e9d 100644 --- a/src/treetensornetworks/opsum_to_ttn.jl +++ b/src/treetensornetworks/opsum_to_ttn.jl @@ -66,6 +66,7 @@ function ordered_incident_edges(g::AbstractGraph,v,ordered_edges) edge_positions=map(e -> findfirst(isequal(e),ordered_edges),edges) return edges[sortperm(edge_positions)] end + # # Tree adaptations of functionalities in ITensors.jl/src/physics/autompo/opsum_to_mpo.jl # @@ -147,123 +148,120 @@ function ttn_svd( #ToDo: precompute span of each term and store # compute span of each term spans=Dict{eltype(os),Vector{vertextype_sites}}() - @time "span" begin for term in os spans[term]=span(term,sites) end - end # build compressed finite state machine representation - @time "first" begin - for v in vs - # for every vertex, find all edges that contain this vertex - edges=ordered_incident_edges(sites,v,es) + for v in vs + # for every vertex, find all edges that contain this vertex + edges=ordered_incident_edges(sites,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] + + which_incident_edge=map_vertices_to_incident_edges(sites,v,edges) + factor_to_edge(t)=which_incident_edge[ITensors.site(t)] + # 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 - # 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] + v in spans[term] || continue + factors=ITensors.terms(term) + + # filter out factor that acts on current vertex + onsite = filter(t -> (ITensors.site(t) == v), factors) + not_onsite_factors=setdiff(factors,onsite) + + # filter out factors that come in from the direction of the incoming edge + incoming=filter(t->factor_to_edge(t)==edge_in, not_onsite_factors) + + # also store all non-incoming factors in standard order, used for channel merging + not_incoming = filter(t-> (ITensors.site(t) == v) || factor_to_edge(t)!=edge_in,factors) - which_incident_edge=map_vertices_to_incident_edges(sites,v,edges) - factor_to_edge(t)=which_incident_edge[ITensors.site(t)] - # 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 every outgoing edge, filter out factors that go out along that edge + outgoing = Dict( + e => filter(t->factor_to_edge(t)==e,not_onsite_factors) + for e in edges_out + ) - for term in os - # loop over OpSum and pick out terms that act on current vertex - - v in spans[term] || continue - factors=ITensors.terms(term) - - # filter out factor that acts on current vertex - onsite = filter(t -> (ITensors.site(t) == v), factors) - not_onsite_factors=setdiff(factors,onsite) - - # filter out factors that come in from the direction of the incoming edge - incoming=filter(t->factor_to_edge(t)==edge_in, not_onsite_factors) - - # also store all non-incoming factors in standard order, used for channel merging - not_incoming = filter(t-> (ITensors.site(t) == v) || factor_to_edge(t)!=edge_in,factors) - - # for every outgoing edge, filter out factors that go out along that edge - outgoing = Dict( - e => filter(t->factor_to_edge(t)==e,not_onsite_factors) - for e in edges_out + # compute QNs + incoming_qn = calc_qn(incoming) + not_incoming_qn = calc_qn(not_incoming) + 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{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 + if !isempty(incoming) + # get the correct map from edge=>QN to term and channel + # this checks if term exists on edge=>QN ( otherwise insert it) and returns it's index + 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 = ITensorMPS.posInLink!(cinmap, incoming) + bond_col = ITensorMPS.posInLink!(coutmap, not_incoming) # get incoming channel + bond_coef = convert(coefficient_type, ITensors.coefficient(term)) + q_inbond_coefs = get!( + inbond_coefs[edge_in], incoming_qn, ITensorMPS.MatElem{coefficient_type}[] ) - - # compute QNs - incoming_qn = calc_qn(incoming) - not_incoming_qn = calc_qn(not_incoming) - 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{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 - if !isempty(incoming) - # get the correct map from edge=>QN to term and channel - # this checks if term exists on edge=>QN ( otherwise insert it) and returns it's index - 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 = ITensorMPS.posInLink!(cinmap, incoming) - bond_col = ITensorMPS.posInLink!(coutmap, not_incoming) # get incoming channel - bond_coef = convert(coefficient_type, ITensors.coefficient(term)) - q_inbond_coefs = get!( - inbond_coefs[edge_in], incoming_qn, ITensorMPS.MatElem{coefficient_type}[] - ) - push!(q_inbond_coefs, ITensorMPS.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 - coutmap = get!( - outmaps, edges[dout] => outgoing_qns[edges[dout]], Dict{Vector{Op},Int}() - ) - T_inds[dout] = ITensorMPS.posInLink!(coutmap, outgoing[edges[dout]]) # add outgoing channel - T_qns[dout] = outgoing_qns[edges[dout]] - end - # if term starts at this site, add its coefficient as a site factor - site_coef = one(coefficient_type) - if (isnothing(dim_in) || T_inds[dim_in] == -1) && - ITensors.argument(term) ∉ site_coef_done - site_coef = ITensors.coefficient(term) - site_coef = convert(coefficient_type, site_coef) # required since ITensors.coefficient seems to return ComplexF64 even if coefficient_type is determined to be real - 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 + push!(q_inbond_coefs, ITensorMPS.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 + coutmap = get!( + outmaps, edges[dout] => outgoing_qns[edges[dout]], Dict{Vector{Op},Int}() + ) + T_inds[dout] = ITensorMPS.posInLink!(coutmap, outgoing[edges[dout]]) # add outgoing channel + T_qns[dout] = outgoing_qns[edges[dout]] + end + # if term starts at this site, add its coefficient as a site factor + site_coef = one(coefficient_type) + if (isnothing(dim_in) || T_inds[dim_in] == -1) && + ITensors.argument(term) ∉ site_coef_done + site_coef = ITensors.coefficient(term) + site_coef = convert(coefficient_type, site_coef) # required since ITensors.coefficient seems to return ComplexF64 even if coefficient_type is determined to be real + 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 - # save indices and value of symbolic tensor entry - el = QNArrElem(T_qns, T_inds, site_coef * Prod(onsite)) - - push!(tempTTN[v], el) end + # save indices and value of symbolic tensor entry + el = QNArrElem(T_qns, T_inds, site_coef * Prod(onsite)) - ITensorMPS.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 = ITensorMPS.toMatrix(mat) - U, S, V = svd(M) - P = S .^ 2 - truncate!(P; maxdim, cutoff, mindim) - tdim = length(P) - nc = size(M, 2) - Vs[edges[dim_in]][q] = Matrix{coefficient_type}(V[1:nc, 1:tdim]) - end + push!(tempTTN[v], el) + end + + ITensorMPS.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 = ITensorMPS.toMatrix(mat) + U, S, V = svd(M) + P = S .^ 2 + truncate!(P; maxdim, cutoff, mindim) + tdim = length(P) + nc = size(M, 2) + Vs[edges[dim_in]][q] = Matrix{coefficient_type}(V[1:nc, 1:tdim]) end end end + # compress this tempTTN representation into dense form link_space = Dict{edgetype_sites,Index}() @@ -554,10 +552,10 @@ function ttn( is_leaf(sites, root_vertex) || error("Tree root must be a leaf vertex.") os = deepcopy(os) - @time os = sorteachterm(os, sites, root_vertex) - @time os = ITensorMPS.sortmergeterms(os) # not exported + os = sorteachterm(os, sites, root_vertex) + os = ITensorMPS.sortmergeterms(os) # not exported if algorithm == "svd" - @time T = ttn_svd(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 From 58503dc4fd24d8d1417fcb2abda5a6ad95fee53b Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 12 Apr 2024 14:51:01 -0400 Subject: [PATCH 05/13] Format. --- src/treetensornetworks/opsum_to_ttn.jl | 95 +++++++++++++------------- 1 file changed, 47 insertions(+), 48 deletions(-) diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl index 7e893e9d..4bbf9116 100644 --- a/src/treetensornetworks/opsum_to_ttn.jl +++ b/src/treetensornetworks/opsum_to_ttn.jl @@ -13,14 +13,12 @@ using StaticArrays: MVector # Utility methods # - - # 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) nterms == 1 && return [ITensors.site(t[1])] - for i in 1:nterms, j in i+1:nterms + for i in 1:nterms, j in (i + 1):nterms path = vertex_path(g, ITensors.site(t[i]), ITensors.site(t[j])) spn = union(spn, path) end @@ -34,36 +32,36 @@ end # map all vertices `w` of `g` except for `v` to the incident edge of `v` # which lies in edge_path(g,w,v) -function map_vertices_to_incident_edges(g::AbstractGraph,v,es) - _g=deepcopy(ITensorNetworks.underlying_graph(g)) - NamedGraphs.rem_vertex!(_g,v) - subgraphs=NamedGraphs.connected_components(_g) - - vs=vertextype(g)[] +function map_vertices_to_incident_edges(g::AbstractGraph, v, es) + _g = deepcopy(ITensorNetworks.underlying_graph(g)) + NamedGraphs.rem_vertex!(_g, v) + subgraphs = NamedGraphs.connected_components(_g) + + vs = vertextype(g)[] for e in es - push!(vs,only(setdiff([src(e),dst(e)],[v,]))) + push!(vs, only(setdiff([src(e), dst(e)], [v]))) end - _vs=vertextype(g)[] - _es=edgetype(g)[] - - for (e,v) in zip(es,vs) + _vs = vertextype(g)[] + _es = edgetype(g)[] + + for (e, v) in zip(es, vs) for i in eachindex(subgraphs) if v in subgraphs[i] - append!(_vs,subgraphs[i]) - append!(_es,fill(e,length(subgraphs[i]))) - deleteat!(subgraphs,i) + append!(_vs, subgraphs[i]) + append!(_es, fill(e, length(subgraphs[i]))) + deleteat!(subgraphs, i) break end end end @assert isempty(subgraphs) - return Dict(zip(_vs,_es)) + return Dict(zip(_vs, _es)) end -function ordered_incident_edges(g::AbstractGraph,v,ordered_edges) - edges=ITensorNetworks.incident_edges(g,v) - edges=map(e -> e in ordered_edges ? e : reverse(e), edges) - edge_positions=map(e -> findfirst(isequal(e),ordered_edges),edges) +function ordered_incident_edges(g::AbstractGraph, v, ordered_edges) + edges = ITensorNetworks.incident_edges(g, v) + edges = map(e -> e in ordered_edges ? e : reverse(e), edges) + edge_positions = map(e -> findfirst(isequal(e), ordered_edges), edges) return edges[sortperm(edge_positions)] end @@ -100,9 +98,9 @@ function ttn_svd( thishasqns = any(v -> hasqns(sites[v]), vertices(sites)) # traverse tree outwards from root vertex - vs = _default_vertex_ordering(sites,root_vertex) + vs = _default_vertex_ordering(sites, root_vertex) # ToDo: Add check in ttn_svd that the ordering matches that of find_index_in_tree, which is used in sorteachterm #fermion-sign! - es = _default_edge_ordering(sites,root_vertex) # store edges in fixed ordering relative to root + es = _default_edge_ordering(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 @@ -147,49 +145,50 @@ function ttn_svd( tempTTN = Dict(v => QNArrElem{Scaled{coefficient_type,Prod{Op}},degrees[v]}[] for v in vs) #ToDo: precompute span of each term and store # compute span of each term - spans=Dict{eltype(os),Vector{vertextype_sites}}() + spans = Dict{eltype(os),Vector{vertextype_sites}}() for term in os - spans[term]=span(term,sites) + spans[term] = span(term, sites) end # build compressed finite state machine representation for v in vs # for every vertex, find all edges that contain this vertex - edges=ordered_incident_edges(sites,v,es) - + edges = ordered_incident_edges(sites, 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] - - which_incident_edge=map_vertices_to_incident_edges(sites,v,edges) - factor_to_edge(t)=which_incident_edge[ITensors.site(t)] + + which_incident_edge = map_vertices_to_incident_edges(sites, v, edges) + factor_to_edge(t) = which_incident_edge[ITensors.site(t)] # 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 - + v in spans[term] || continue - factors=ITensors.terms(term) + factors = ITensors.terms(term) # filter out factor that acts on current vertex onsite = filter(t -> (ITensors.site(t) == v), factors) - not_onsite_factors=setdiff(factors,onsite) + not_onsite_factors = setdiff(factors, onsite) # filter out factors that come in from the direction of the incoming edge - incoming=filter(t->factor_to_edge(t)==edge_in, not_onsite_factors) - + incoming = filter(t -> factor_to_edge(t) == edge_in, not_onsite_factors) + # also store all non-incoming factors in standard order, used for channel merging - not_incoming = filter(t-> (ITensors.site(t) == v) || factor_to_edge(t)!=edge_in,factors) - + not_incoming = filter( + t -> (ITensors.site(t) == v) || factor_to_edge(t) != edge_in, factors + ) + # for every outgoing edge, filter out factors that go out along that edge outgoing = Dict( - e => filter(t->factor_to_edge(t)==e,not_onsite_factors) - for e in edges_out + e => filter(t -> factor_to_edge(t) == e, not_onsite_factors) for e in edges_out ) - + # compute QNs incoming_qn = calc_qn(incoming) not_incoming_qn = calc_qn(not_incoming) @@ -261,7 +260,7 @@ function ttn_svd( end end end - + # compress this tempTTN representation into dense form link_space = Dict{edgetype_sites,Index}() @@ -288,7 +287,7 @@ function ttn_svd( for v in vs # redo the whole thing like before # ToDo: use neighborhood instead of going through all edges, see above - edges=ordered_incident_edges(sites,v,es) + edges = ordered_incident_edges(sites, v, es) dim_in = findfirst(e -> dst(e) == v, edges) dims_out = findall(e -> src(e) == v, edges) # slice isometries at this vertex @@ -452,11 +451,11 @@ function computeSiteProd(sites::IndsNetwork{V,<:Index}, ops::Prod{Op})::ITensor return T end -function _default_vertex_ordering(g::AbstractGraph,root_vertex) +function _default_vertex_ordering(g::AbstractGraph, root_vertex) return reverse(post_order_dfs_vertices(g, root_vertex)) end -function _default_edge_ordering(g::AbstractGraph,root_vertex) +function _default_edge_ordering(g::AbstractGraph, root_vertex) return reverse(reverse.(post_order_dfs_edges(g, root_vertex))) end @@ -467,8 +466,8 @@ function sorteachterm(os::OpSum, sites::IndsNetwork{V,<:Index}, root_vertex::V) os = copy(os) # linear ordering of vertices in tree graph relative to chosen root, chosen outward from root - ordering = _default_vertex_ordering(sites,root_vertex) - site_positions=Dict(zip(ordering,1:length(ordering))) + ordering = _default_vertex_ordering(sites, root_vertex) + site_positions = Dict(zip(ordering, 1:length(ordering))) findpos(op::Op) = site_positions[ITensors.site(op)] isless_site(o1::Op, o2::Op) = findpos(o1) < findpos(o2) N = nv(sites) From b97da4ee143a6c760f13ae8f78b2e55c7c65a674 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 12 Apr 2024 14:57:54 -0400 Subject: [PATCH 06/13] Fix some unnecessary namespace qualifiers and change copy to deepcopy. --- 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 4bbf9116..ebcdae59 100644 --- a/src/treetensornetworks/opsum_to_ttn.jl +++ b/src/treetensornetworks/opsum_to_ttn.jl @@ -33,9 +33,9 @@ end # map all vertices `w` of `g` except for `v` to the incident edge of `v` # which lies in edge_path(g,w,v) function map_vertices_to_incident_edges(g::AbstractGraph, v, es) - _g = deepcopy(ITensorNetworks.underlying_graph(g)) - NamedGraphs.rem_vertex!(_g, v) - subgraphs = NamedGraphs.connected_components(_g) + _g = copy(underlying_graph(g)) + rem_vertex!(_g, v) + subgraphs = connected_components(_g) vs = vertextype(g)[] for e in es @@ -59,7 +59,7 @@ function map_vertices_to_incident_edges(g::AbstractGraph, v, es) end function ordered_incident_edges(g::AbstractGraph, v, ordered_edges) - edges = ITensorNetworks.incident_edges(g, v) + edges = incident_edges(g, v) edges = map(e -> e in ordered_edges ? e : reverse(e), edges) edge_positions = map(e -> findfirst(isequal(e), ordered_edges), edges) return edges[sortperm(edge_positions)] From 7147c6e393167bed057ab4e3a1961dfa0fed4a93 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 12 Apr 2024 15:41:12 -0400 Subject: [PATCH 07/13] Convert some Vectors to Sets. --- src/treetensornetworks/opsum_to_ttn.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl index ebcdae59..f73d7ca9 100644 --- a/src/treetensornetworks/opsum_to_ttn.jl +++ b/src/treetensornetworks/opsum_to_ttn.jl @@ -15,11 +15,11 @@ using StaticArrays: MVector # determine 'support' of product operator on tree graph function span(t::Scaled{C,Prod{Op}}, g::AbstractGraph) where {C} - spn = eltype(g)[] + spn = Set{eltype(g)}() nterms = length(t) - nterms == 1 && return [ITensors.site(t[1])] + nterms == 1 && return Set([ITensors.site(t[1])]) for i in 1:nterms, j in (i + 1):nterms - path = vertex_path(g, ITensors.site(t[i]), ITensors.site(t[j])) + path = Set(vertex_path(g, ITensors.site(t[i]), ITensors.site(t[j]))) spn = union(spn, path) end return spn @@ -35,7 +35,7 @@ end function map_vertices_to_incident_edges(g::AbstractGraph, v, es) _g = copy(underlying_graph(g)) rem_vertex!(_g, v) - subgraphs = connected_components(_g) + subgraphs = Set.(connected_components(_g)) vs = vertextype(g)[] for e in es @@ -145,7 +145,7 @@ function ttn_svd( tempTTN = Dict(v => QNArrElem{Scaled{coefficient_type,Prod{Op}},degrees[v]}[] for v in vs) #ToDo: precompute span of each term and store # compute span of each term - spans = Dict{eltype(os),Vector{vertextype_sites}}() + spans = Dict{eltype(os),Set{vertextype_sites}}() for term in os spans[term] = span(term, sites) end From 9f7b2aba0ecd3f22a8a6a9f49586577770769c77 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 12 Apr 2024 18:04:30 -0400 Subject: [PATCH 08/13] Addressed comments. --- src/treetensornetworks/opsum_to_ttn.jl | 43 ++++++++++++++------------ 1 file changed, 23 insertions(+), 20 deletions(-) diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl index f73d7ca9..68c1e916 100644 --- a/src/treetensornetworks/opsum_to_ttn.jl +++ b/src/treetensornetworks/opsum_to_ttn.jl @@ -20,7 +20,7 @@ function span(t::Scaled{C,Prod{Op}}, g::AbstractGraph) where {C} nterms == 1 && return Set([ITensors.site(t[1])]) for i in 1:nterms, j in (i + 1):nterms path = Set(vertex_path(g, ITensors.site(t[i]), ITensors.site(t[j]))) - spn = union(spn, path) + spn = union!(spn, path) end return spn end @@ -30,21 +30,28 @@ function crosses_vertex(t::Scaled{C,Prod{Op}}, g::AbstractGraph, v) where {C} return v ∈ span(t, g) end -# map all vertices `w` of `g` except for `v` to the incident edge of `v` +function align_edges(reference_edges, edges) + return intersect(reference_edges, Iterators.flatten((edges, reverse.(edges)))) +end + +# return a dict from vertices `w` of `g`, except for `v`, to the incident edge of `v` # which lies in edge_path(g,w,v) -function map_vertices_to_incident_edges(g::AbstractGraph, v, es) +function vertices_to_incident_edges_dict(g::AbstractGraph, v, incident_edges) + #split graph into subtrees by removing vertex v _g = copy(underlying_graph(g)) rem_vertex!(_g, v) subgraphs = Set.(connected_components(_g)) + #for each incident edge, store the vertex that's not `v` vs = vertextype(g)[] - for e in es + for e in incident_edges push!(vs, only(setdiff([src(e), dst(e)], [v]))) end + + #return a Dictionary from vertices to incident_edges _vs = vertextype(g)[] _es = edgetype(g)[] - - for (e, v) in zip(es, vs) + for (e, v) in zip(incident_edges, vs) for i in eachindex(subgraphs) if v in subgraphs[i] append!(_vs, subgraphs[i]) @@ -58,13 +65,6 @@ function map_vertices_to_incident_edges(g::AbstractGraph, v, es) return Dict(zip(_vs, _es)) end -function ordered_incident_edges(g::AbstractGraph, v, ordered_edges) - edges = incident_edges(g, v) - edges = map(e -> e in ordered_edges ? e : reverse(e), edges) - edge_positions = map(e -> findfirst(isequal(e), ordered_edges), edges) - return edges[sortperm(edge_positions)] -end - # # Tree adaptations of functionalities in ITensors.jl/src/physics/autompo/opsum_to_mpo.jl # @@ -152,7 +152,7 @@ function ttn_svd( # build compressed finite state machine representation for v in vs # for every vertex, find all edges that contain this vertex - edges = ordered_incident_edges(sites, v, es) + edges = align_edges(es, incident_edges(sites, v)) # use the corresponding ordering as index order for tensor elements at this site dim_in = findfirst(e -> dst(e) == v, edges) @@ -160,8 +160,7 @@ function ttn_svd( dims_out = findall(e -> src(e) == v, edges) edges_out = edges[dims_out] - which_incident_edge = map_vertices_to_incident_edges(sites, v, edges) - factor_to_edge(t) = which_incident_edge[ITensors.site(t)] + which_incident_edge = vertices_to_incident_edges_dict(sites, v, edges) # 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) @@ -177,16 +176,20 @@ function ttn_svd( not_onsite_factors = setdiff(factors, onsite) # filter out factors that come in from the direction of the incoming edge - incoming = filter(t -> factor_to_edge(t) == edge_in, not_onsite_factors) + incoming = filter( + t -> which_incident_edge[ITensors.site(t)] == edge_in, not_onsite_factors + ) # also store all non-incoming factors in standard order, used for channel merging not_incoming = filter( - t -> (ITensors.site(t) == v) || factor_to_edge(t) != edge_in, factors + t -> (ITensors.site(t) == v) || which_incident_edge[ITensors.site(t)] != edge_in, + factors, ) # for every outgoing edge, filter out factors that go out along that edge outgoing = Dict( - e => filter(t -> factor_to_edge(t) == e, not_onsite_factors) for e in edges_out + e => filter(t -> which_incident_edge[ITensors.site(t)] == e, not_onsite_factors) for + e in edges_out ) # compute QNs @@ -287,7 +290,7 @@ function ttn_svd( for v in vs # redo the whole thing like before # ToDo: use neighborhood instead of going through all edges, see above - edges = ordered_incident_edges(sites, v, es) + edges = align_edges(es, incident_edges(sites, v)) dim_in = findfirst(e -> dst(e) == v, edges) dims_out = findall(e -> src(e) == v, edges) # slice isometries at this vertex From 283e9bd4c24f2af8ddc062c0cb16d89b70314639 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Fri, 12 Apr 2024 18:07:28 -0400 Subject: [PATCH 09/13] Change argumetn order for align_edges. --- 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 68c1e916..f7bad79e 100644 --- a/src/treetensornetworks/opsum_to_ttn.jl +++ b/src/treetensornetworks/opsum_to_ttn.jl @@ -30,7 +30,7 @@ function crosses_vertex(t::Scaled{C,Prod{Op}}, g::AbstractGraph, v) where {C} return v ∈ span(t, g) end -function align_edges(reference_edges, edges) +function align_edges(edges, reference_edges) return intersect(reference_edges, Iterators.flatten((edges, reverse.(edges)))) end @@ -152,7 +152,7 @@ function ttn_svd( # build compressed finite state machine representation for v in vs # for every vertex, find all edges that contain this vertex - edges = align_edges(es, incident_edges(sites, v)) + edges = align_edges(incident_edges(sites, v), es) # use the corresponding ordering as index order for tensor elements at this site dim_in = findfirst(e -> dst(e) == v, edges) @@ -290,7 +290,7 @@ function ttn_svd( for v in vs # redo the whole thing like before # ToDo: use neighborhood instead of going through all edges, see above - edges = align_edges(es, incident_edges(sites, v)) + edges = align_edges(incident_edges(sites, v), es) dim_in = findfirst(e -> dst(e) == v, edges) dims_out = findall(e -> src(e) == v, edges) # slice isometries at this vertex From cba04e71a681b892f4841708f63b47bfd1bced47 Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Mon, 15 Apr 2024 16:05:01 -0400 Subject: [PATCH 10/13] Remove use of span, move some newly introduced functions into ttn_svd. --- src/treetensornetworks/opsum_to_ttn.jl | 93 ++++++++++---------------- 1 file changed, 37 insertions(+), 56 deletions(-) diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl index f7bad79e..7c9dd383 100644 --- a/src/treetensornetworks/opsum_to_ttn.jl +++ b/src/treetensornetworks/opsum_to_ttn.jl @@ -6,65 +6,30 @@ using ITensors.NDTensors: Block, maxdim, nblocks, nnzblocks using ITensors.Ops: Op, OpSum using NamedGraphs: degrees, is_leaf, vertex_path using StaticArrays: MVector - +using DataGraphs: AbstractDataGraph +using NamedGraphs: boundary_edges # convert ITensors.OpSum to TreeTensorNetwork # # Utility methods # -# determine 'support' of product operator on tree graph -function span(t::Scaled{C,Prod{Op}}, g::AbstractGraph) where {C} - spn = Set{eltype(g)}() - nterms = length(t) - nterms == 1 && return Set([ITensors.site(t[1])]) - for i in 1:nterms, j in (i + 1):nterms - path = Set(vertex_path(g, ITensors.site(t[i]), ITensors.site(t[j]))) - spn = union!(spn, path) - end - return spn +function align_edges(edges, reference_edges) + return intersect(Iterators.flatten(zip(edges, reverse.(edges))), reference_edges) 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) +function align_and_reorder_edges(edges, reference_edges) + return intersect(reference_edges, align_edges(edges, reference_edges)) end -function align_edges(edges, reference_edges) - return intersect(reference_edges, Iterators.flatten((edges, reverse.(edges)))) -end - -# return a dict from vertices `w` of `g`, except for `v`, to the incident edge of `v` -# which lies in edge_path(g,w,v) -function vertices_to_incident_edges_dict(g::AbstractGraph, v, incident_edges) - #split graph into subtrees by removing vertex v - _g = copy(underlying_graph(g)) +function split_at_vertex(g::AbstractGraph, v) + _g = copy(g) rem_vertex!(_g, v) - subgraphs = Set.(connected_components(_g)) - - #for each incident edge, store the vertex that's not `v` - vs = vertextype(g)[] - for e in incident_edges - push!(vs, only(setdiff([src(e), dst(e)], [v]))) - end - - #return a Dictionary from vertices to incident_edges - _vs = vertextype(g)[] - _es = edgetype(g)[] - for (e, v) in zip(incident_edges, vs) - for i in eachindex(subgraphs) - if v in subgraphs[i] - append!(_vs, subgraphs[i]) - append!(_es, fill(e, length(subgraphs[i]))) - deleteat!(subgraphs, i) - break - end - end - end - @assert isempty(subgraphs) - return Dict(zip(_vs, _es)) + return Set.(connected_components(_g)) end +split_at_vertex(g::AbstractDataGraph, v) = split_at_vertex(underlying_graph(g), v) + # # Tree adaptations of functionalities in ITensors.jl/src/physics/autompo/opsum_to_mpo.jl # @@ -143,16 +108,11 @@ function ttn_svd( 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) - #ToDo: precompute span of each term and store - # compute span of each term - spans = Dict{eltype(os),Set{vertextype_sites}}() - for term in os - spans[term] = span(term, sites) - end + # build compressed finite state machine representation for v in vs # for every vertex, find all edges that contain this vertex - edges = align_edges(incident_edges(sites, v), es) + edges = align_and_reorder_edges(incident_edges(sites, v), es) # use the corresponding ordering as index order for tensor elements at this site dim_in = findfirst(e -> dst(e) == v, edges) @@ -160,7 +120,19 @@ function ttn_svd( dims_out = findall(e -> src(e) == v, edges) edges_out = edges[dims_out] - which_incident_edge = vertices_to_incident_edges_dict(sites, v, edges) + # for every site w except v, determine the incident edge to v that lies + # in the edge_path(w,v) + subgraphs = split_at_vertex(sites, v) + _boundary_edges = align_edges( + [only(boundary_edges(underlying_graph(sites), subgraph)) for subgraph in subgraphs], + edges, + ) + which_incident_edge = Dict( + Iterators.flatten([ + subgraphs[i] .=> ((_boundary_edges[i]),) for i in eachindex(subgraphs) + ]), + ) + # 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) @@ -168,8 +140,17 @@ function ttn_svd( for term in os # loop over OpSum and pick out terms that act on current vertex - v in spans[term] || continue factors = ITensors.terms(term) + if v in ITensors.site.(factors) + crosses_vertex = true + else + crosses_vertex = + !isone( + length(Set([which_incident_edge[site] for site in ITensors.site.(factors)])) + ) + end + #if term doesn't cross vertex, skip it + crosses_vertex || continue # filter out factor that acts on current vertex onsite = filter(t -> (ITensors.site(t) == v), factors) @@ -290,7 +271,7 @@ function ttn_svd( for v in vs # redo the whole thing like before # ToDo: use neighborhood instead of going through all edges, see above - edges = align_edges(incident_edges(sites, v), es) + edges = align_and_reorder_edges(incident_edges(sites, v), es) dim_in = findfirst(e -> dst(e) == v, edges) dims_out = findall(e -> src(e) == v, edges) # slice isometries at this vertex From 2b734634e2ad20ceb4ac11a283fb98887d50322d Mon Sep 17 00:00:00 2001 From: Benedikt Kloss Date: Tue, 16 Apr 2024 11:01:12 -0400 Subject: [PATCH 11/13] Bump compat entry for DataGraphs and remove dispatch on AbstractDataGraphs for connected_components. --- Project.toml | 4 ++-- src/treetensornetworks/opsum_to_ttn.jl | 5 +---- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/Project.toml b/Project.toml index 2ce895f0..f9893795 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ITensorNetworks" uuid = "2919e153-833c-4bdc-8836-1ea460a35fc7" authors = ["Matthew Fishman and contributors"] -version = "0.7" +version = "0.7.0" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" @@ -45,7 +45,7 @@ ITensorNetworksEinExprsExt = "EinExprs" AbstractTrees = "0.4.4" Combinatorics = "1" Compat = "3, 4" -DataGraphs = "0.1.7" +DataGraphs = "v0.1.13" DataStructures = "0.18" Dictionaries = "0.4" Distributions = "0.25.86" diff --git a/src/treetensornetworks/opsum_to_ttn.jl b/src/treetensornetworks/opsum_to_ttn.jl index 7c9dd383..7fd22851 100644 --- a/src/treetensornetworks/opsum_to_ttn.jl +++ b/src/treetensornetworks/opsum_to_ttn.jl @@ -6,7 +6,6 @@ using ITensors.NDTensors: Block, maxdim, nblocks, nnzblocks using ITensors.Ops: Op, OpSum using NamedGraphs: degrees, is_leaf, vertex_path using StaticArrays: MVector -using DataGraphs: AbstractDataGraph using NamedGraphs: boundary_edges # convert ITensors.OpSum to TreeTensorNetwork @@ -28,8 +27,6 @@ function split_at_vertex(g::AbstractGraph, v) return Set.(connected_components(_g)) end -split_at_vertex(g::AbstractDataGraph, v) = split_at_vertex(underlying_graph(g), v) - # # Tree adaptations of functionalities in ITensors.jl/src/physics/autompo/opsum_to_mpo.jl # @@ -37,7 +34,7 @@ split_at_vertex(g::AbstractDataGraph, v) = split_at_vertex(underlying_graph(g), """ ttn_svd(os::OpSum, sites::IndsNetwork, root_vertex, kwargs...) -Construct a dense TreeTensorNetwork from a symbolic OpSum representation of a +Construct a TreeTensorNetwork from a symbolic OpSum representation of a Hamiltonian, compressing shared interaction channels. """ function ttn_svd(os::OpSum, sites::IndsNetwork, root_vertex; kwargs...) From 05ffb5d9c6e1749f4e51d1426d03439b7d5cc251 Mon Sep 17 00:00:00 2001 From: b-kloss Date: Tue, 16 Apr 2024 11:14:23 -0400 Subject: [PATCH 12/13] Update Project.toml Co-authored-by: Matt Fishman --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ae38f99a..ff911a1e 100644 --- a/Project.toml +++ b/Project.toml @@ -46,7 +46,7 @@ ITensorNetworksEinExprsExt = "EinExprs" AbstractTrees = "0.4.4" Combinatorics = "1" Compat = "3, 4" -DataGraphs = "v0.1.13" +DataGraphs = "0.1.13" DataStructures = "0.18" Dictionaries = "0.4" Distributions = "0.25.86" From 618e298788ab8a8dfb4d8c70e7d5e696361d88e9 Mon Sep 17 00:00:00 2001 From: b-kloss Date: Tue, 16 Apr 2024 11:14:29 -0400 Subject: [PATCH 13/13] Update src/treetensornetworks/opsum_to_ttn.jl Co-authored-by: Matt Fishman --- 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 7fd22851..a3923dd8 100644 --- a/src/treetensornetworks/opsum_to_ttn.jl +++ b/src/treetensornetworks/opsum_to_ttn.jl @@ -22,9 +22,9 @@ function align_and_reorder_edges(edges, reference_edges) end function split_at_vertex(g::AbstractGraph, v) - _g = copy(g) - rem_vertex!(_g, v) - return Set.(connected_components(_g)) + g = copy(g) + rem_vertex!(g, v) + return Set.(connected_components(g)) end #