Kd/use factor graphs (#56)

* some attempts to use pepsrow

* test 4x4 grid

* correction to case 2

* correction

* test on the number of vertices

* correction

* correction

* corrections

* some change

* some change

* correction

* dicts

* tests on tensor

* tests on grit formation

* a picture

* corrections

* corrections

* make a test more clear

* add test for pathological instance

* move reshape_row to utils.jl

* add test on reshape_row

* passing tests

* changes

* add test file

* tests that does not work

* correction

* correction

* some tests on chimera

* some corrections

* solved problem with lattice

* corrections

Co-authored-by: Bartłomiej Gardas <bartek.gardas@gmail.com>
Co-authored-by: annamariadziubyna <73058800+annamariadziubyna@users.noreply.github.com>

benchmarks/tests_grid_from_file.jl

@@ -65,14 +65,23 @@ for k in 1:examples
     energies_given = data["energies"][k,:,:]
 
     g = M2graph(Mat_of_interactions, -1)
+    si = Int(sqrt(size(Mat_of_interactions, 1)))
 
     ################ exact method ###################
-
+    fg = factor_graph(
+        g,
+        2,
+        energy=energy,
+        spectrum=brute_force,
+    )
+
+    peps = PepsNetwork(si, si, fg, β, :NW)
+    println("size of peps = ", peps.size)
     print("peps time  = ")
 
 
     number = number_of_states + more_states_for_peps
-    @time spins, objective = solve(g, number ; β = T(β), threshold = 0., δH = δH)
+    @time spins, objective = solve(g, peps, number ; β = T(β), threshold = 0., δH = δH)
 
     for i in 1:number_of_states
 
@@ -88,7 +97,7 @@ for k in 1:examples
     print("approx peps  ")
 
     number = number_of_states + more_states_for_peps
-    @time spins_approx, objective_approx = solve(g, number; β = T(β), χ = χ, threshold = 1e-12, δH = δH)
+    @time spins_approx, objective_approx = solve(g, peps, number; β = T(β), χ = χ, threshold = 1e-12, δH = δH)
 
     for i in 1:number_of_states
 
@@ -102,8 +111,43 @@ for k in 1:examples
     @test objective ≈ objective_approx atol = 1.e-7
 
     print("peps larger T")
+    s1 = ceil(Int, si/2)
+    s2 = floor(Int, si/2)
+    println(s1)
+
+    if s1 == s2
+        rule = square_lattice((s1, 2, s1, 2, 1))
+    else
+        D1 = Dict{Any,Any}(1 => 1, 2 => 1, 6 => 1, 7 => 1)
+        D2 = Dict{Any,Any}(3 => 2, 4 => 2, 8 => 2, 9 => 2)
+        D3 = Dict{Any,Any}(5 => 3, 10 => 3)
+        D4 = Dict{Any,Any}(11 => 4, 12 => 4, 16 => 4, 17 => 4)
+        D5 = Dict{Any,Any}(13 => 5, 14 => 5, 18 => 5, 19 => 5)
+        D6 = Dict{Any,Any}(15 => 6, 20 => 6)
+        D7 = Dict{Any,Any}(21 => 7, 22 => 7)
+        D8 = Dict{Any,Any}(23 => 8, 24 => 8)
+        D9 = Dict{Any,Any}(25 => 9)
+
+        rule = merge(D1, D2, D3, D4, D5, D6, D7, D8, D9)
+    end
+
+    update_cells!(
+      g,
+      rule = rule,
+    )
+
+    fg = factor_graph(
+        g,
+        16,
+        energy=energy,
+        spectrum=brute_force,
+    )
+
+    peps = PepsNetwork(s1, s1, fg, β, :NW)
+    println("size of peps = ", peps.size)
+
     number = number_of_states + more_states_for_peps
-    @time spins_larger_nodes, objective_larger_nodes = solve(g, number; node_size = (2,2), β = T(β), χ = χ, threshold = 1e-12, δH = δH)
+    @time spins_larger_nodes, objective_larger_nodes = solve(g, peps, number; node_size = (2,2), β = T(β), χ = χ, threshold = 1e-12, δH = δH)
 
     for i in 1:number_of_states
 
@@ -117,8 +161,16 @@ for k in 1:examples
     @test objective ≈ objective_larger_nodes atol = 1.e-7
 
     print("peps larger T, limited spectrum")
+    fg = factor_graph(
+        g,
+        15,
+        energy=energy,
+        spectrum=brute_force,
+    )
+
+    peps = PepsNetwork(s1, s1, fg, β, :NW)
     number = number_of_states + more_states_for_peps
-    @time spins_spec, objective_spec = solve(g, number; node_size = (2,2), β = T(β), χ = χ, threshold = 1e-12, spectrum_cutoff = 15, δH = δH)
+    @time spins_spec, objective_spec = solve(g, peps, number; node_size = (2,2), β = T(β), χ = χ, threshold = 1e-12, spectrum_cutoff = 15, δH = δH)
 
     for i in 1:minimum([number_of_states, 60])
         @test energy(spins_spec[i], g) ≈ energies_given[i]

benchmarks/tests_on_chimera.jl

@@ -54,7 +54,7 @@ s = ArgParseSettings("description")
     arg_type = Int
     help = "chimera node size in columns"
     "--spectrum_cutoff", "-u"
-    default = 1000
+    default = 256
     arg_type = Int
     help = "size of the lower spectrum"
     "--deltaH", "-d"
@@ -74,15 +74,37 @@ problem_size = parse_args(s)["size"]
 χ = parse_args(s)["chi"]
 si = parse_args(s)["size"]
 δH = parse_args(s)["deltaH"]
+spectrum_cutoff = parse_args(s)["spectrum_cutoff"]
+node_size = (parse_args(s)["node_size1"], parse_args(s)["node_size2"])
+s1 = Int(sqrt(si/8))
+
+n = ceil(Int, s1/node_size[1])
+m = ceil(Int, s1/node_size[2])
 
 ig = ising_graph(fi, si, 1)
 
+ig1 = ising_graph(fi, si, 1)
+update_cells!(
+    ig1,
+    rule = square_lattice((m, node_size[1], n, node_size[2], 8)),
+  )
+
+fg = factor_graph(
+      ig1,
+      spectrum_cutoff,
+      energy=energy,
+      spectrum=brute_force,
+  )
+
+peps = PepsNetwork(m, n, fg, β, :NW)
+
+println("size od peps = ", peps.size)
+
 n_sols = parse_args(s)["n_sols"]
-node_size = (parse_args(s)["node_size1"], parse_args(s)["node_size2"])
+
 println(node_size)
-spectrum_cutoff = parse_args(s)["spectrum_cutoff"]
 
-@time spins, objective = solve(ig, n_sols; β=β, χ = χ, threshold = 1e-8, node_size = node_size, spectrum_cutoff = spectrum_cutoff, δH=δH)
+@time spins, objective = solve(ig, peps, n_sols; β=β, χ = χ, threshold = 1e-8, node_size = node_size, spectrum_cutoff = spectrum_cutoff, δH=δH)
 
 energies = [energy(s, ig) for s in spins]
 println("energies from peps")

src/lattice.jl

@@ -1,7 +1,7 @@
 export square_lattice
 
  function square_lattice(size::NTuple{5, Int})  
-    m, um, n, un, t = size  
+    m, um, n, un, t = size
     new = LinearIndices((1:n, 1:m))
     old = LinearIndices((1:t, 1:un, 1:n, 1:um, 1:m))
 
@@ -12,7 +12,7 @@ export square_lattice
     rule
 end
 
-function square_lattice(size::NTuple{3, Int})  
-    m, n, t = size  
+function square_lattice(size::NTuple{3, Int})
+    m, n, t = size
     square_lattice((m, 1, n, 1, t))
-end
+end

src/peps_no_types.jl

@@ -352,12 +352,13 @@ function merge_dX(partial_s::Vector{Partial_sol{T}}, dX_inds::Vector{Int}, δH::
 end
 
 
-function solve(g::MetaGraph, no_sols::Int = 2; node_size::Tuple{Int, Int} = (1,1),
+function solve(g::MetaGraph, peps::PepsNetwork, no_sols::Int = 2; node_size::Tuple{Int, Int} = (1,1),
                                                β::T, χ::Int = 2^prod(node_size),
                                                threshold::Float64 = 0.,
                                                spectrum_cutoff::Int = 1000,
                                                δH::Float64 = 0.) where T <: Real
 
+
     gg = graph4peps(g, node_size, spectrum_cutoff = spectrum_cutoff)
 
     grid = props(gg)[:grid]

src/utils.jl

@@ -52,7 +52,7 @@ local_basis(d::Int) = union(-1, 1:d-1)
 local_basis(ψ::AbstractMPS, i::Int) = local_basis(size(ψ[i], 2))
 
 function proj(state, dims::Union{Vector, NTuple})
-    P = Matrix{Float64}[] 
+    P = Matrix{Float64}[]
     for (σ, r) ∈ zip(state, dims)
         v = zeros(r)
         v[idx(σ)...] = 1.
@@ -126,7 +126,7 @@ function Base.LinearIndices(m::Int, n::Int, origin::Symbol=:NW)
 
     if origin ∈ (:NW, :NE, :SE, :SW)
         i_max, j_max = m, n
-    else 
+    else
         i_max, j_max = n, m
     end
 
@@ -217,4 +217,4 @@ function rank_vec(ig::MetaGraph)
     rank = get_prop(ig, :rank)
     L = get_prop(ig, :L)
     Int[get(rank, i, 1) for i=1:L]
-end
+end

test/PEPS.jl

@@ -114,7 +114,6 @@ instance = "$(@__DIR__)/instances/$(L)_001.txt"
 
 ig = ising_graph(instance, L)
 
-#states = collect.(all_states(get_prop(ig, :rank)))
 states = collect.(all_states(rank_vec(ig)))
 ρ = exp.(-β .* energy.(states, Ref(ig)))
 Z = sum(ρ)

test/base.jl

@@ -18,13 +18,13 @@ T = ComplexF64
     @test rank(ψ) == Tuple(fill(d, 1:sites))
     @test bond_dimension(ψ) ≈ D
 
-    ϕ = copy(ψ) 
+    ϕ = copy(ψ)
     @test ϕ == ψ
     @test ϕ ≈ ψ
 
     show(ψ)
 
-    dims = (3, 2, 5, 4) 
+    dims = (3, 2, 5, 4)
     @info "Veryfing ψ of arbitrary rank" dims
 
     ψ = randn(MPS{T}, D, dims)
@@ -39,7 +39,7 @@ T = ComplexF64
     @test rank(ψ) == dims
     @test bond_dimension(ψ) ≈ D
 
-    ϕ = copy(ψ) 
+    ϕ = copy(ψ)
     @test ϕ == ψ
     @test ϕ ≈ ψ
 
@@ -57,40 +57,20 @@ end
     @test eltype(O) == ComplexF64
 
     P = copy(O)
-    @test P == O 
-    @test P ≈ O 
+    @test P == O
+    @test P ≈ O
 end
 
-@testset "Reshaping (row-wise)" begin
-    vec = Vector(1:6)
-
-    A = reshape_row(vec, (2, 3))
-    B = [1 2 3; 4 5 6]
-
-    @test A ≈ B
-end 
-
-@testset "Basic vector to tensor reshaping" begin
-    dims = (2, 3, 4, 5)
-    states = [randn(T, d) for d ∈ dims] 
-    vec = kron(states...)
-
-    ψ = tensor(MPS(states))
-    ϕ = reshape_row(vec, dims)
-
-    @test ψ ≈ ϕ
-end 
-
 @testset "MPS from tensor" begin
     ϵ = 1E-14
 
-    dims = (2,3,4,3,5) 
+    dims = (2,3,4,3,5)
     sites = length(dims)
-    A = randn(T, dims) 
+    A = randn(T, dims)
 
     @test sqrt(sum(abs.(A) .^ 2)) ≈ norm(A)
 
-    @test ndims(A) == sites 
+    @test ndims(A) == sites
     @test size(A) == dims
 
     ψ = MPS(A, :right)
@@ -106,18 +86,18 @@ end
     @test rank(ψ) == size(AA)
     @test norm(AA) ≈ 1
     @test size(AA) == size(A)
- 
+
     vA = vec(A)
     nA = norm(vA)
-    @test abs(1 - abs(dot(vec(AA), vA ./ nA))) < ϵ 
+    @test abs(1 - abs(dot(vec(AA), vA ./ nA))) < ϵ
     #@test AA ≈ A ./ norm(A) # this is true "module phase"
-    
+
     B = randn(T, dims...)
     ϕ = MPS(B, :left)
 
     @test norm(ϕ) ≈ 1
     @test_nowarn verify_bonds(ϕ)
-    @test_nowarn verify_physical_dims(ϕ, dims)    
+    @test_nowarn verify_physical_dims(ϕ, dims)
     @test is_left_normalized(ϕ)
     show(ϕ)
 
@@ -126,16 +106,16 @@ end
     @test rank(ϕ) == size(BB)
     @test norm(BB) ≈ 1
     @test sqrt(sum(abs.(B) .^ 2)) ≈ norm(B)
-    
+
     vB = vec(B)
     nB = norm(vB)
-    @test abs(1 - abs(dot(vec(BB), vB ./ nB))) < ϵ     
+    @test abs(1 - abs(dot(vec(BB), vB ./ nB))) < ϵ
     #@test BB ≈ B ./ norm(B) # this is true "module phase"
-    
+
     χ = MPS(A, :left)
 
     @test norm(χ) ≈ 1
-    @test abs(1 - abs(dot(ψ, χ))) < ϵ 
+    @test abs(1 - abs(dot(ψ, χ))) < ϵ
 end
 
-end
+end

test/graph.jl

@@ -124,11 +124,16 @@ update_cells!(
 
 fg = factor_graph(
     ig,
-
     energy=energy,
     spectrum=full_spectrum,
 )
 
+for v ∈ vertices(fg)
+   cl = get_prop(fg, v, :cluster)
+   nodes = [e for e in keys(cl.vertices)]
+   @test sort(nodes) == cells[v]
+end
+
 
 for (bd, e) in zip(bond_dimensions, edges(fg))
    pl, en, pr = get_prop(fg, e, :split)

test/mps_tests.jl

@@ -91,7 +91,16 @@ end
 
     spins_exact, _ = solve_mps(g, 10; β=β, β_step=1, χ=12, threshold = 0.)
 
-    spins_peps, _ = solve(g, 10; β = β, χ = 2, threshold = 1e-12, δH = 0.1)
+    fg = factor_graph(
+        g,
+        2,
+        energy=energy,
+        spectrum=brute_force,
+    )
+
+    peps = PepsNetwork(4, 4, fg, β, :NW)
+
+    spins_peps, _ = solve(g, peps, 10; β = β, χ = 2, threshold = 1e-12, δH = 0.1)
 
     for k in 1:10
         #testing exact