some corrections in merging solutions

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))
-δ = 0.9
+δH = 0.9
 # this is axiliary function for npz write
 
 function vecvec2matrix(v::Vector{Vector{Int}})
@@ -101,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])
 

benchmarks/tests_grid_from_file.jl

@@ -31,7 +31,7 @@ println("examples = ", examples)
 
 
 β = 3.
-δ = 0.
+δH = 0.
 
 number_of_states = 10
 
@@ -72,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
 
@@ -88,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
 
@@ -103,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
 
@@ -118,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]

benchmarks/tests_on_chimera.jl

@@ -23,7 +23,6 @@ function vecvec2matrix(v::Vector{Vector{Int}})
     end
     M
 end
-δ = 0.1
 
 s = ArgParseSettings("description")
   @add_arg_table! s begin
@@ -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,6 +73,7 @@ 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)
 
@@ -78,7 +82,7 @@ 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_no_types.jl

@@ -292,18 +292,71 @@ function conditional_probabs(gg::MetaGraph, ps::Partial_sol{T}, j::Int, lower_mp
     probs_unnormed./sum(probs_unnormed)
 end
 
+"""
+    function dX_inds(grid::Matrix{Int}, j::Int; has_diagonals::Bool = false)
+
+Returns vector{Int} indexing of the boundary region (dX) given a grid.
+id has diagonals, diagonal bounds on the grid are taken into account
+"""
 
-function indices_on_boundary(grid::Matrix{Int}, j::Int)
-    smax = j-1
-    smin = maximum([1, j - size(grid, 2)])
-    collect(smin: smax)
+function dX_inds(grid::Matrix{Int}, j::Int; has_diagonals::Bool = false)
+    last = j-1
+    s = size(grid, 2)
+    first = maximum([1, j - s])
+    if (has_diagonals & (j%s != 1))
+        first = maximum([1, j - s - 1])
+    end
+    return collect(first: last)
 end
 
+"""
+    function merge_dX(partial_s::Vector{Partial_sol{T}}, dX_inds::Vector{Int}, δH::Float64) where T <:Real
+
+Return a vector of Partial_sol{T}, with merged boundaries.
+
+Merging rule is such that the retion of the objective function of the merged item
+to the maximal is lower than δH
+"""
+
+function merge_dX(partial_s::Vector{Partial_sol{T}}, dX_inds::Vector{Int}, δH::Float64) where T <:Real
+    if (length(partial_s) > 1) & (δH != .0)
+        leave = [true for _ in partial_s]
+
+        dXes = [ps.spins[dX_inds] for ps in partial_s]
+
+        unique_dXes = unique(dXes)
+        if dXes != unique_dXes
+            dXcount = countmap(dXes)
+            for dX in unique_dXes
+                if dXcount[dX] > 1
+                    i = findall(k -> k == dX, dXes)
+                    objectives = [partial_s[j].objective for j in i]
+
+                    objectives = objectives./maximum(objectives)
+                    for ind in i[objectives .< δH]
+                        leave[ind] = false
+                    end
+                end
+            end
+            no_reduced = count(.!(leave))
+            # this is just for testing
+            if no_reduced > 0
+                j = length(partial_s[1].spins)
+                k = length(partial_s)
+                println(no_reduced, " out of $k partial solutions deleted at j = $j")
+            end
+            return partial_s[leave]
+        end
+    end
+    partial_s
+end
+
+
 function solve(g::MetaGraph, no_sols::Int = 2; node_size::Tuple{Int, Int} = (1,1),
                                                β::T, χ::Int = 2^prod(node_size),
                                                threshold::Float64 = 0.,
                                                spectrum_cutoff::Int = 1000,
-                                               δ::Float64 = 0.1) where T <: Real
+                                               δH::Float64 = 0.) where T <: Real
 
     gg = graph4peps(g, node_size, spectrum_cutoff = spectrum_cutoff)
 
@@ -319,11 +372,11 @@ function solve(g::MetaGraph, no_sols::Int = 2; node_size::Tuple{Int, Int} = (1,1
 
         for j in grid[row,:]
 
-            dX = indices_on_boundary(grid, j)
+            dX = dX_inds(grid, j)
 
             partial_s_temp = Partial_sol{T}[]
 
-            partial_s = merge_boundaries(partial_s, dX, δ)
+            partial_s = merge_dX(partial_s, dX, δH)
             for ps in partial_s
 
                 objectives = conditional_probabs(gg, ps, j, lower_mps, vec_of_T)
@@ -378,33 +431,6 @@ function return_solutions(partial_s::Vector{Partial_sol{T}}, ns:: MetaGraph)  wh
     return spins, objective
 end
 
-function merge_boundaries(partial_s::Vector{Partial_sol{T}}, dX::Vector{Int}, δ) where T <:Real
-    if (length(partial_s) > 1) & (δ != .0)
-        leave = [true for _ in partial_s]
-
-        boundaries = [ps.spins[dX] for ps in partial_s]
-
-        unique_bondaries = unique(boundaries)
-        if boundaries != unique_bondaries
-            b = countmap(boundaries)
-            for el in unique_bondaries
-                if b[el] > 1
-                    i = findall(x -> x == el, boundaries)
-                    objectives = [partial_s[j].objective for j in i]
-
-                    objectives = objectives./maximum(objectives)
-                    for ind in i[objectives .< δ]
-                        leave[ind] = false
-                        println(leave)
-                    end
-                end
-            end
-            return partial_s[leave]
-        end
-    end
-    partial_s
-end
-
 """
     select_best_solutions(partial_s_temp::Vector{Partial_sol{T}}, no_sols::Int) where T <:Real
 

test/mps_tests.jl

@@ -91,7 +91,7 @@ end
 
     spins_exact, _ = solve_mps(g, 10; β=β, β_step=1, χ=12, threshold = 0.)
 
-    spins_peps, _ = solve(g, 10; β = β, χ = 2, threshold = 1e-12)
+    spins_peps, _ = solve(g, 10; β = β, χ = 2, threshold = 1e-12, δH = 0.1)
 
     for k in 1:10
         #testing exact