Add definitions of delta and Ising networks (#18)

Author: mtfishman

Project.toml

@@ -1,13 +1,14 @@
 name = "ITensorNetworksNext"
 uuid = "302f2e75-49f0-4526-aef7-d8ba550cb06c"
 authors = ["ITensor developers <[email protected]> and contributors"]
-version = "0.1.12"
+version = "0.1.13"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 BackendSelection = "680c2d7c-f67a-4cc9-ae9c-da132b1447a5"
 DataGraphs = "b5a273c3-7e6c-41f6-98bd-8d7f1525a36a"
+DiagonalArrays = "74fd4be6-21e2-4f6f-823a-4360d37c7a77"
 Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -31,6 +32,7 @@ AbstractTrees = "0.4.5"
 Adapt = "4.3"
 BackendSelection = "0.1.6"
 DataGraphs = "0.2.7"
+DiagonalArrays = "0.3.23"
 Dictionaries = "0.4.5"
 Graphs = "1.13.1"
 LinearAlgebra = "1.10"

docs/src/reference.md

@@ -1,5 +1,5 @@
 # Reference
 
 ```@autodocs
-Modules = [ITensorNetworksNext]
+Modules = [ITensorNetworksNext, ITensorNetworksNext.TensorNetworkGenerators]
 ```

src/ITensorNetworksNext.jl

@@ -3,6 +3,7 @@ module ITensorNetworksNext
 include("lazynameddimsarrays.jl")
 include("abstracttensornetwork.jl")
 include("tensornetwork.jl")
+include("TensorNetworkGenerators/TensorNetworkGenerators.jl")
 include("contract_network.jl")
 include("abstract_problem.jl")
 include("iterators.jl")

src/TensorNetworkGenerators/TensorNetworkGenerators.jl

@@ -0,0 +1,8 @@
+module TensorNetworkGenerators
+
+export delta_network, ising_network
+
+include("delta_network.jl")
+include("ising_network.jl")
+
+end

src/TensorNetworkGenerators/delta_network.jl

@@ -0,0 +1,21 @@
+using DiagonalArrays: δ
+using Graphs: AbstractGraph
+using ..ITensorNetworksNext: TensorNetwork
+using NamedGraphs.GraphsExtensions: incident_edges
+
+"""
+    delta_network(f, elt::Type = Float64, g::AbstractGraph)
+
+Construct a TensorNetwork on the graph `g` with element type `elt` that has delta tensors
+on each vertex. Link dimensions are defined using the function `f(e)` that should take an
+edge `e` as an input and should output the link index on that edge.
+"""
+function delta_network(f, elt::Type, g::AbstractGraph)
+    return tn = TensorNetwork(g) do v
+        is = Tuple(f.(incident_edges(g, v)))
+        return δ(elt, is)
+    end
+end
+function delta_network(f, g::AbstractGraph)
+    return delta_network(f, Float64, g)
+end

src/TensorNetworkGenerators/ising_network.jl

@@ -0,0 +1,47 @@
+using DiagonalArrays: DiagonalArray
+using Graphs: degree, dst, edges, src
+using LinearAlgebra: Diagonal, eigen
+using NamedDimsArrays: apply, dename, inds, operator, randname
+
+function sqrt_ising_bond(β; h1 = zero(typeof(β)), h2 = zero(typeof(β)))
+    f11 = exp(β * (1 + h1 + h2))
+    f12 = exp(β * (-1 + h1 - h2))
+    f21 = exp(β * (-1 - h1 + h2))
+    f22 = exp(β * (1 - h1 - h2))
+    m² = eltype(β)[f11 f12; f21 f22]
+    d², v = eigen(m²)
+    d = sqrt.(d²)
+    return v * Diagonal(d) * inv(v)
+end
+
+"""
+    ising_network(f, β::Number, g::AbstractGraph)
+
+Construct a TensorNetwork on the graph `g` with inverse temperature `β` that has Ising
+partition function tensors on each vertex. Link dimensions are defined using the function
+`f(e)` that should take an edge `e` as an input and should output the link index on that
+edge.
+"""
+function ising_network(
+        f, β::Number, g::AbstractGraph; h::Number = zero(eltype(β)), sz_vertices = []
+    )
+    elt = typeof(β)
+    l̃ = Dict(e => randname(f(e)) for e in edges(g))
+    f̃(e) = get(() -> l̃[reverse(e)], l̃, e)
+    tn = delta_network(f̃, elt, g)
+    for v in sz_vertices
+        a = DiagonalArray(elt[1, -1], dename.(inds(tn[v])))
+        tn[v] = a[inds(tn[v])...]
+    end
+    for e in edges(tn)
+        v1 = src(e)
+        v2 = dst(e)
+        deg1 = degree(tn, v1)
+        deg2 = degree(tn, v2)
+        m = sqrt_ising_bond(β; h1 = h / deg1, h2 = h / deg2)
+        t = operator(m, (f̃(e),), (f(e),))
+        tn[v1] = apply(t, tn[v1])
+        tn[v2] = apply(t, tn[v2])
+    end
+    return tn
+end

src/abstracttensornetwork.jl

@@ -1,37 +1,16 @@
 using Adapt: Adapt, adapt, adapt_structure
 using BackendSelection: @Algorithm_str, Algorithm
-using DataGraphs:
-    DataGraphs,
-    AbstractDataGraph,
-    edge_data,
-    underlying_graph,
-    underlying_graph_type,
-    vertex_data
+using DataGraphs: DataGraphs, AbstractDataGraph, edge_data, underlying_graph,
+    underlying_graph_type, vertex_data
 using Dictionaries: Dictionary
-using Graphs:
-    Graphs,
-    AbstractEdge,
-    AbstractGraph,
-    Graph,
-    add_edge!,
-    add_vertex!,
-    bfs_tree,
-    center,
-    dst,
-    edges,
-    edgetype,
-    ne,
-    neighbors,
-    nv,
-    rem_edge!,
-    src,
-    vertices
+using Graphs: Graphs, AbstractEdge, AbstractGraph, Graph, add_edge!, add_vertex!,
+    bfs_tree, center, dst, edges, edgetype, ne, neighbors, nv, rem_edge!, src, vertices
 using LinearAlgebra: LinearAlgebra, factorize
 using MacroTools: @capture
 using NamedDimsArrays: dimnames, inds
 using NamedGraphs: NamedGraphs, NamedGraph, not_implemented, steiner_tree
-using NamedGraphs.GraphsExtensions:
-    ⊔, directed_graph, incident_edges, rem_edges!, rename_vertices, vertextype
+using NamedGraphs.GraphsExtensions: ⊔, directed_graph, incident_edges, rem_edges!,
+    rename_vertices, vertextype
 using SplitApplyCombine: flatten
 
 abstract type AbstractTensorNetwork{V, VD} <: AbstractDataGraph{V, VD, Nothing} end
@@ -202,7 +181,6 @@ end
 # Fix the edges of the TensorNetwork `tn` to match
 # the tensor connectivity at vertex `v`.
 function fix_edges!(tn::AbstractTensorNetwork, v)
-    rem_incident_edges!(tn, v)
     rem_edges!(tn, incident_edges(tn, v))
     add_missing_edges!(tn, v)
     return tn

test/Project.toml

@@ -1,12 +1,14 @@
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+DiagonalArrays = "74fd4be6-21e2-4f6f-823a-4360d37c7a77"
 Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7"
 ITensorNetworksNext = "302f2e75-49f0-4526-aef7-d8ba550cb06c"
 NamedDimsArrays = "60cbd0c0-df58-4cb7-918c-6f5607b73fde"
 NamedGraphs = "678767b0-92e7-4007-89e4-4527a8725b19"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
@@ -17,15 +19,17 @@ WrappedUnions = "325db55a-9c6c-5b90-b1a2-ec87e7a38c44"
 [compat]
 AbstractTrees = "0.4.5"
 Aqua = "0.8.14"
+DiagonalArrays = "0.3.23"
 Dictionaries = "0.4.5"
 Graphs = "1.13.1"
 ITensorBase = "0.3"
 ITensorNetworksNext = "0.1.1"
 NamedDimsArrays = "0.8"
 NamedGraphs = "0.6.8, 0.7"
+QuadGK = "2.11.2"
 SafeTestsets = "0.1"
 Suppressor = "0.2.8"
-TermInterface = "2"
 TensorOperations = "5.3.1"
+TermInterface = "2"
 Test = "1.10"
 WrappedUnions = "0.3"

test/test_tensornetworkgenerators.jl

@@ -0,0 +1,62 @@
+using DiagonalArrays: δ
+using Graphs: edges, ne, nv, vertices
+using ITensorBase: Index
+using ITensorNetworksNext: contract_network
+using ITensorNetworksNext.TensorNetworkGenerators: delta_network, ising_network
+using NamedDimsArrays: inds
+using NamedGraphs.GraphsExtensions: arranged_edges, incident_edges
+using NamedGraphs.NamedGraphGenerators: named_grid
+using Test: @test, @testset
+
+module TestUtils
+    using QuadGK: quadgk
+    # Exact critical inverse temperature for 2D square lattice Ising model.
+    βc() = 0.5 * log(1 + √2)
+    # Exact infinite volume free energy density for 2D square lattice Ising model.
+    function ising_free_energy_density(β::Real)
+        κ = 2sinh(2β) / cosh(2β)^2
+        integrand(θ) = log(0.5 * (1 + sqrt(abs(1 - (κ * sin(θ))^2))))
+        integral, _ = quadgk(integrand, 0, π)
+        return (-log(2cosh(2β)) - (1 / (2π)) * integral) / β
+    end
+end
+
+@testset "TensorNetworkGenerators" begin
+    @testset "Delta Network" begin
+        dims = (3, 3)
+        g = named_grid(dims)
+        ldict = Dict(e => Index(2) for e in edges(g))
+        l(e) = get(() -> ldict[reverse(e)], ldict, e)
+        tn = delta_network(l, g)
+        @test nv(tn) == 9
+        @test ne(tn) == ne(g)
+        @test issetequal(vertices(tn), vertices(g))
+        @test issetequal(arranged_edges(tn), arranged_edges(g))
+        for v in vertices(tn)
+            is = l.(incident_edges(g, v))
+            @test tn[v] == δ(Tuple(is))
+        end
+    end
+    @testset "Ising Network" begin
+        dims = (4, 4)
+        β = TestUtils.βc()
+        g = named_grid(dims; periodic = true)
+        ldict = Dict(e => Index(2) for e in edges(g))
+        l(e) = get(() -> ldict[reverse(e)], ldict, e)
+        tn = ising_network(l, β, g)
+        @test nv(tn) == 16
+        @test ne(tn) == ne(g)
+        @test issetequal(vertices(tn), vertices(g))
+        @test issetequal(arranged_edges(tn), arranged_edges(g))
+        for v in vertices(tn)
+            is = l.(incident_edges(g, v))
+            @test issetequal(is, inds(tn[v]))
+            @test tn[v] ≠ δ(Tuple(is))
+        end
+        # TODO: Use eager contraction sequence finding.
+        z = contract_network(tn; alg = "exact")[]
+        f = -log(z) / (β * nv(g))
+        f_inf = TestUtils.ising_free_energy_density(β)
+        @test f ≈ f_inf rtol = 1.0e-1
+    end
+end