Kd graphic

Project.toml

@@ -7,13 +7,19 @@ version = "0.0.1"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b"
 LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 LowRankApprox = "898213cb-b102-5a47-900c-97e73b919f73"
 MetaGraphs = "626554b9-1ddb-594c-aa3c-2596fe9399a5"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TensorCast = "02d47bb6-7ce6-556a-be16-bb1710789e2b"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [compat]
 MetaGraphs = "0.6"
@@ -25,6 +31,7 @@ LowRankApprox = "0.4"
 Requires = "1.1"
 TensorCast = "0.3"
 TensorOperations = "3.0.1"
+StatsBase = "0.33"
 julia = "1.5"
 
 [extras]

benchmarks/low_energy_approx_chimera.jl

@@ -14,7 +14,7 @@ using Test
 import SpinGlassPEPS: solve, solve_mps, M2graph, energy, binary2spins, ising_graph
 
 disable_logging(LogLevel(0))
-
+δH = 0.9
 # this is axiliary function for npz write
 
 function vecvec2matrix(v::Vector{Vector{Int}})
@@ -28,7 +28,7 @@ end
 
 s = ArgParseSettings("description")
   @add_arg_table! s begin
-    "--file", "-f"
+    "--file", "-i"
     arg_type = String
     help = "the file name"
     "--size", "-s"
@@ -72,7 +72,7 @@ problem_size = parse_args(s)["size"]
 χ = parse_args(s)["chi"]
 si = parse_args(s)["size"]
 
-ig = ising_graph(fi, si, 1, 1)
+ig = ising_graph(fi, si, 1)
 
 n_sols = parse_args(s)["n_sols"]
 node_size = (parse_args(s)["node_size1"], parse_args(s)["node_size2"])
@@ -88,9 +88,8 @@ energy_ref = energy(ground_spins, ig)
 
 spectrum_cutoff = parse_args(s)["spectrum_cutoff"]
 
-ses = collect(spectrum_cutoff:-10:1)
-step = 9
-step = 100
+ses = collect(spectrum_cutoff:-10:40)
+step = 10
 n_s = collect(n_sols:-step:1)
 
 delta_e = ones(length(ses), length(n_s))
@@ -102,7 +101,7 @@ function proceed()
     i = 1
     for s in ses
 
-      @time spins, _ = solve(ig, n_sol; β=β, χ = χ, threshold = 1e-8, node_size = node_size, spectrum_cutoff = s)
+      @time spins, _ = solve(ig, n_sol; β=β, χ = χ, threshold = 1e-8, node_size = node_size, spectrum_cutoff = s, δH=δH)
 
       en = minimum([energy(s, ig) for s in spins])
 

test/tests_grid_from_file.jl → benchmarks/tests_grid_from_file.jl

@@ -31,6 +31,7 @@ println("examples = ", examples)
 
 
 β = 3.
+δH = 0.
 
 number_of_states = 10
 
@@ -71,7 +72,7 @@ for k in 1:examples
 
 
     number = number_of_states + more_states_for_peps
-    @time spins, objective = solve(g, number ; β = T(β), threshold = 0.)
+    @time spins, objective = solve(g, number ; β = T(β), threshold = 0., δH = δH)
 
     for i in 1:number_of_states
 
@@ -87,7 +88,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)
+    @time spins_approx, objective_approx = solve(g, number; β = T(β), χ = χ, threshold = 1e-12, δH = δH)
 
     for i in 1:number_of_states
 
@@ -102,7 +103,7 @@ for k in 1:examples
 
     print("peps larger T")
     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)
+    @time spins_larger_nodes, objective_larger_nodes = solve(g, number; node_size = (2,2), β = T(β), χ = χ, threshold = 1e-12, δH = δH)
 
     for i in 1:number_of_states
 
@@ -117,7 +118,7 @@ for k in 1:examples
 
     print("peps larger T, limited spectrum")
     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)
+    @time spins_spec, objective_spec = solve(g, 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]

test/tests_on_chimera.jl → benchmarks/tests_on_chimera.jl

@@ -24,10 +24,9 @@ function vecvec2matrix(v::Vector{Vector{Int}})
     M
 end
 
-
 s = ArgParseSettings("description")
   @add_arg_table! s begin
-    "--file", "-f"
+    "--file", "-i"
     arg_type = String
     help = "the file name"
     "--size", "-s"
@@ -58,6 +57,10 @@ s = ArgParseSettings("description")
     default = 1000
     arg_type = Int
     help = "size of the lower spectrum"
+    "--deltaH", "-d"
+    default = 0.1
+    arg_type = Float64
+    help = "merge parameter on merging dX, the threshold on ratios of objectives"
   end
 
 fi = parse_args(s)["file"]
@@ -70,15 +73,16 @@ problem_size = parse_args(s)["size"]
 β = parse_args(s)["beta"]
 χ = parse_args(s)["chi"]
 si = parse_args(s)["size"]
+δH = parse_args(s)["deltaH"]
 
-ig = ising_graph(fi, si, 1, 1)
+ig = ising_graph(fi, si, 1)
 
 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)
+@time spins, objective = solve(ig, 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/PEPS.jl

@@ -1,65 +1,166 @@
-export PepsTensor
-
-mutable struct PepsTensor
-    nbrs::Dict{String, Int}
-    left::AbstractArray
-    right::AbstractArray
-    up::AbstractArray
-    down::AbstractArray
-    loc::AbstractArray
-    tensor::AbstractArray
-
-    function PepsTensor(fg::MetaDiGraph, v::Int)
-        pc = new()
-        pc.nbrs = Dict()
-        pc.loc = get_prop(fg, v, :local_exp)
-
-        outgoing = outneighbors(fg, v)
-        incoming = inneighbors(fg, v)
-
-        for u ∈ outgoing 
-            e = SimpleEdge(v, u)
-            if get_prop(fg, e, :orientation) == "horizontal"
-                pc.right = first(get_prop(fg, e, :decomposition))
-                push!(pc.nbrs, "h_out" => u)
-            else
-                pc.down = first(get_prop(fg, e, :decomposition))
-                push!(pc.nbrs, "v_out" => u)
-            end 
-        end
+export NetworkGraph, PepsNetwork
+export generate_tensor, MPO
+
+mutable struct NetworkGraph
+    factor_graph::MetaDiGraph
+    nbrs::Dict
+    β::Number
+
+    function NetworkGraph(factor_graph::MetaDiGraph, nbrs::Dict, β::Number)
+        ng = new(factor_graph, nbrs, β)
 
-        for u ∈ incoming
-            e = SimpleEdge(u, v)
-            if get_prop(fg, e, :orientation) == "horizontal"
-                pc.left = last(get_prop(fg, e, :decomposition))
-                push!(pc.nbrs, "h_in" => u)
-            else
-                pc.up = last(get_prop(fg, e, :decomposition))
-                push!(pc.nbrs, "v_in" => u)
-            end 
+        count = 0
+        for v ∈ vertices(ng.factor_graph), w ∈ ng.nbrs[v]
+            if has_edge(ng.factor_graph, v, w) count += 1 end
         end
-       
-        # open boundary conditions
-        if !isdefined(pc, :left)
-            pc.left = ones(1, size(pc.right, 1))
+
+        mc = ne(ng.factor_graph)
+        if count < mc
+            error("Error: $(count) < $(mc)")
         end
+        ng
+    end
+end
+
+function generate_tensor(ng::NetworkGraph, v::Int)
+    fg = ng.factor_graph
+    loc_exp = exp.(-ng.β .* get_prop(fg, v, :loc_en))
+
+    dim = []
+    @cast tensor[_, i] := loc_exp[i]
 
-        if !isdefined(pc, :right)
-            pc.right = ones(size(pc.left, 2), 1)
+    for w ∈ ng.nbrs[v]
+        if has_edge(fg, w, v)
+            _, _, pv = get_prop(fg, w, v, :split)
+            pv = pv'
+        elseif has_edge(fg, v, w)
+            pv, _, _ = get_prop(fg, v, w, :split)
+        else
+            pv = ones(length(loc_exp), 1)
         end
 
-        if !isdefined(pc, :up)
-            pc.up = ones(1, size(pc.down, 1))
+        @cast tensor[(c, γ), σ] |= tensor[c, σ] * pv[σ, γ]
+        push!(dim, size(pv, 2))
+    end
+
+    reshape(tensor, dim..., :)
+end
+
+function _generate_tensor(ng::NetworkGraph, v::Int)
+    fg = ng.factor_graph
+    loc_exp = exp.(-ng.β .* get_prop(fg, v, :loc_en))
+
+    projs = Dict()
+    for (i, w) ∈ enumerate(ng.nbrs[v])
+        if has_edge(fg, w, v)
+            _, _, pv = get_prop(fg, w, v, :split)
+            pv = pv'
+        elseif has_edge(fg, v, w)
+            pv, _, _ = get_prop(fg, v, w, :split)
+        else
+            pv = ones(length(loc_exp), 1)
         end
+        push!(projs, i => pv)
+    end
+
+    L, U, R, D = projs[1], projs[2], projs[3], projs[4]
+    @cast tensor[l, u, r, d, σ] |= L[σ, l] * U[σ, u] * R[σ, r] * D[σ, d] * loc_exp[σ]
 
-        if !isdefined(pc, :down)
-            pc.down = ones(size(pc.up, 2), 1)
+    tensor
+end
+
+function generate_tensor(ng::NetworkGraph, v::Int, w::Int)
+    fg = ng.factor_graph
+    if has_edge(fg, w, v)
+        _, e, _ = get_prop(fg, w, v, :split)
+        tensor = exp.(-ng.β .* e')
+    elseif has_edge(fg, v, w)
+        _, e, _ = get_prop(fg, v, w, :split)
+        tensor = exp.(-ng.β .* e)
+    else
+        tensor = ones(1, 1)
+    end
+    tensor
+end
+
+mutable struct PepsNetwork
+    size::NTuple{2, Int}
+    map::Dict
+    network_graph::NetworkGraph
+    origin::Symbol
+    i_max::Int
+    j_max::Int
+
+    function PepsNetwork(m::Int, n::Int, fg::MetaDiGraph, β::Number, origin::Symbol)
+        pn = new((m, n))
+        pn.map, pn.i_max, pn.j_max = LinearIndices(m, n, origin)
+
+        nbrs = Dict()
+        for i ∈ 1:pn.i_max, j ∈ 1:pn.j_max
+            push!(nbrs,
+            pn.map[i, j] => (pn.map[i, j-1], pn.map[i-1, j],
+                             pn.map[i, j+1], pn.map[i+1, j]))
         end
+        pn.network_graph = NetworkGraph(fg, nbrs, β)
+        pn
+    end
+end
+
+generate_tensor(pn::PepsNetwork, m::NTuple{2,Int}) = generate_tensor(pn.network_graph, pn.map[m])
+_generate_tensor(pn::PepsNetwork, m::NTuple{2,Int}) = _generate_tensor(pn.network_graph, pn.map[m])
+generate_tensor(pn::PepsNetwork, m::NTuple{2,Int}, n::NTuple{2,Int}) = generate_tensor(pn.network_graph, pn.map[m], pn.map[n])
+
+function MPO(::Type{T}, ψ::PEPSRow) where {T <: Number}
+    n = length(ψ)
+    ϕ = MPO(T, n)
+    for i=1:n
+        A = ψ[i]
+        @reduce B[l, u, r, d] |= sum(σ) A[l, u, r, d, σ]
+        ϕ[i] = B
+    end
+    ϕ
+end
+MPO(ψ::PEPSRow) = MPO(Float64, ψ)
 
-        @cast pc.tensor[l, r, u, d, σ] |= pc.loc[σ] * pc.left[l, σ] * pc.right[σ, r] * pc.up[u, σ] * pc.down[σ, d]
+function PEPSRow(::Type{T}, peps::PepsNetwork, i::Int) where {T <: Number}
+    n = peps.j_max
+    ψ = PEPSRow(T, n)
 
-        pc
+    for j ∈ 1:n
+        ψ[j] = generate_tensor(peps, (i, j))
     end
+
+    for j ∈ 1:n-1
+        ten = generate_tensor(peps, (i, j), (i, j+1))
+        A = ψ[j]
+        @tensor B[l, u, r, d, σ] := A[l, u, r̃, d, σ] * ten[r̃, r]
+         ψ[j] = B
+    end
+    ψ
 end
+PEPSRow(peps::PepsNetwork, i::Int) = PEPSRow(Float64, peps, i)
+
+function MPO(::Type{T}, peps::PepsNetwork, i::Int, k::Int) where {T <: Number}
+    n = peps.j_max
+
+    ψ = MPO(T, n)
+    fg = peps.network_graph.factor_graph
+
+    for j ∈ 1:n
+        v, w = peps.map[i, j], peps.map[k, j]
 
-Base.size(A::PepsTensor) = size(A.tensor)
+        if has_edge(fg, v, w)
+            _, en, _ = get_prop(fg, v, w, :split)
+        elseif has_edge(fg, w, v)
+            _, en, _ = get_prop(fg, w, v, :split)
+        else
+            en = ones(1, 1)
+        end
+
+        #@cast A[_, σ, _, η] := en[σ, η]
+        @cast A[u, _, d, _] := en[u, d]
+        ψ[j] = A
+    end
+    ψ
+end
+MPO(peps::PepsNetwork, i::Int, k::Int) = MPO(Float64, peps, i, k)