first implemntation of mps mpo on the full graph

src/mps_implementation.jl

@@ -1,4 +1,4 @@
-include("peps.jl")
+include("peps_no_types.jl")
 include("notation.jl")
 
 
@@ -195,13 +195,14 @@ function cluster_conncetions(all_is::Vector{Int}, all_js::Vector{Vector{Int}})
     end
 end
 
-function connections_for_mps(grid::Matrix{Int})
+function connections_for_mps(ns::Vector{Node_of_grid})
     all_is = Int[]
     all_js = Vector{Int}[]
     conections = Vector{Int}[]
-    for i in grid
-        left, right, up, down = read_connecting_pairs(grid, i)
-        pairs = [left..., right..., up..., down...]
+
+    for i in [e.i for e in ns]
+
+        pairs = ns[i].all_connections
         pairs_not_accounted = Int[]
         for p in pairs
 
@@ -240,14 +241,12 @@ function construct_mps(qubo::Vector{Qubo_el{T}}, β::T, β_step::Int, l::Int,
 end
 
 
-function solve_mps(qubo::Vector{Qubo_el{T}}, grid::Matrix{Int}, problem_size::Int,
+function solve_mps(qubo::Vector{Qubo_el{T}}, ns::Vector{Node_of_grid},
                 no_sols::Int; β::T, β_step::Int, χ::Int = 0, threshold::T = T(0.)) where T <: AbstractFloat
 
-    physical_dim = 2
-
-    ns = [Node_of_grid(i, grid) for i in 1:maximum(grid)]
 
-    is, js = connections_for_mps(Array(transpose(grid)))
+    problem_size = length(ns)
+    is, js = connections_for_mps(ns)
     all_is, all_js = cluster_conncetions(is,js)
     mps = construct_mps(qubo, β, β_step, problem_size, all_is, all_js, χ, threshold)
 
@@ -258,7 +257,7 @@ function solve_mps(qubo::Vector{Qubo_el{T}}, grid::Matrix{Int}, problem_size::In
         for ps in partial_s
             objectives = compute_probs(mps, ps.spins)
 
-            for l in 1:physical_dim
+            for l in 1:length(objectives)
                 push!(partial_s_temp, add_spin(ps, l, objectives[l]))
             end
         end

src/notation.jl

@@ -19,6 +19,34 @@ function read_connecting_pairs(grid::Matrix{Int}, i::Int)
     return left, right, up, down
 end
 
+struct Bond_with_other_node
+    node::Int
+    spins1::Vector{Int}
+    spins2::Vector{Int}
+end
+
+struct Qubo_el{T<:AbstractFloat}
+    ind::Tuple{Int, Int}
+    coupling::T
+    function(::Type{Qubo_el{T}})(ind::Tuple{Int, Int}, coupling::T1) where {T <: AbstractFloat, T1 <: AbstractFloat}
+        new{T}(ind, T(coupling))
+    end
+    function(::Type{Qubo_el})(ind::Tuple{Int, Int}, coupling::Float64)
+        new{Float64}(ind, coupling)
+    end
+end
+
+function matrix2qubo_vec(M::Matrix{T}) where T <:AbstractFloat
+    s = size(M)
+    q_vec = Qubo_el{T}[]
+    for i in 1:s[1]
+        for j in 1:i
+            push!(q_vec, Qubo_el((i,j), M[i,j]))
+        end
+    end
+    q_vec 
+end
+
 """
     struct Node_of_grid
 
@@ -33,28 +61,43 @@ struct Node_of_grid
     right::Vector{Vector{Int}}
     up::Vector{Vector{Int}}
     down::Vector{Vector{Int}}
+    all_connections::Vector{Vector{Int}}
+    bonds::Vector{Bond_with_other_node}
     function(::Type{Node_of_grid})(i::Int, grid::Matrix{Int})
         s = size(grid)
         intra_struct = Vector{Int}[]
 
         left, right, up, down = read_connecting_pairs(grid, i)
-        new(i, [i], intra_struct, left, right, up, down)
+        all_connections = [left..., right..., up..., down...]
+        #all_connections = [el[2] for el in all_connections]
+        new(i, [i], intra_struct, left, right, up, down, all_connections)
+    end
+    function(::Type{Node_of_grid})(i::Int, qubos::Vector{Qubo_el{T}}) where T <: AbstractFloat
+        x = Vector{Int}[]
+        all_connections = Vector{Int}[]
+        for q in qubos
+            if (q.ind[1] == i && q.ind[2] != i)
+                push!(all_connections, [q.ind...])
+            end
+            if (q.ind[2] == i && q.ind[1] != i)
+                push!(all_connections, [q.ind[2], q.ind[1]])
+            end
+        end
+
+        new(i, [i], x, x, x, x, x, all_connections)
     end
 end
 
-function get_system_size(ns::Vector{Node_of_grid})
-    mapreduce(x -> length(x.spin_inds), +, ns)
+function get_system_size(qubo::Vector{Qubo_el{T}}) where T <: AbstractFloat
+    size = 0
+    for q in qubo
+        size = maximum([size, q.ind[1], q.ind[2]])
+    end
+    size
 end
 
-struct Qubo_el{T<:AbstractFloat}
-    ind::Tuple{Int, Int}
-    coupling::T
-    function(::Type{Qubo_el{T}})(ind::Tuple{Int, Int}, coupling::T1) where {T <: AbstractFloat, T1 <: AbstractFloat}
-        new{T}(ind, T(coupling))
-    end
-    function(::Type{Qubo_el})(ind::Tuple{Int, Int}, coupling::Float64)
-        new{Float64}(ind, coupling)
-    end
+function get_system_size(ns::Vector{Node_of_grid})
+    mapreduce(x -> length(x.spin_inds), +, ns)
 end
 
 

src/tests/mps_tests.jl

@@ -6,7 +6,9 @@ using TensorOperations
     @test b == 16.0*ones(1,1)
 
     M = [1 2; 3 4]
-    i,j = connections_for_mps(Array(transpose(M)))
+    ns = [Node_of_grid(i, Array(transpose(M))) for i in 1:maximum(M)]
+
+    i,j = connections_for_mps(ns)
     @test i == [1,2,3]
     @test j == [[3, 2], [4], [4]]
     a,b = cluster_conncetions(i,j)
@@ -15,17 +17,19 @@ using TensorOperations
     @test b == [[[3, 2]], [[4]], [[4]]]
 
     grid = [1 2 3; 4 5 6; 7 8 9]
+    A = Array(transpose(grid))
 
-    i,j = connections_for_mps(grid)
-    @test i == [1, 4, 7, 2, 5, 8, 3, 6]
-    @test j == [[2, 4], [5, 7], [8], [3, 5], [6, 8], [9], [6], [9]]
+    ns = [Node_of_grid(i, grid) for i in 1:maximum(grid)]
+    i,j = connections_for_mps(ns)
 
-    a,b = cluster_conncetions(i,j)
-    @test a == [[1, 7], [4, 8], [2, 6], [5], [3]]
-    @test b == [[[2, 4], [8]], [[5, 7], [9]], [[3, 5], [9]], [[6, 8]], [[6]]]
+    @test i == [1, 2, 3, 4, 5, 6, 7, 8]
+    @test j == [[2, 4], [3, 5], [6], [5, 7], [6, 8], [9], [8], [9]]
 
+    a,b = cluster_conncetions(i,j)
+    @test a == [[1, 5], [2, 6], [3, 7], [4, 8]]
+    @test b == [[[2, 4], [6, 8]], [[3, 5], [9]], [[6], [8]], [[5, 7], [9]]]
 end
-if true
+
 function make_qubo_x()
     qubo = [(1,1) .5; (1,2) -0.5; (1,4) -1.5; (2,2) -1.; (2,3) -1.5; (2,5) -0.5; (3,3) 2.; (3,6) 1.5]
     qubo = vcat(qubo, [(6,6) .05; (5,6) -0.25; (6,9) -0.52; (5,5) 0.75; (4,5) 0.5; (5,8) 0.5; (4,4) 0.; (4,7) -0.01])
@@ -90,7 +94,9 @@ end
     @test mps[3] == ones(1,1,2)
 
     grid = [1 2 3; 4 5 6; 7 8 9]
-    is, js = connections_for_mps(Array(transpose(grid)))
+    ns = [Node_of_grid(i,grid) for i in 1:maximum(grid)]
+
+    is, js = connections_for_mps(ns)
     all_is, all_js = cluster_conncetions(is,js)
 
     mps = construct_mps(qubo, β, 2, 9, all_is, all_js, 4, 0.)
@@ -164,7 +170,9 @@ end
 
     grid = [1 2 3; 4 5 6; 7 8 9]
 
-    spins, _ = solve_mps(train_qubo, grid, 9, 2; β=2., β_step=2, χ=2, threshold = 1e-14)
+    ns = [Node_of_grid(i, train_qubo) for i in 1:get_system_size(train_qubo)]
+
+    spins, _ = solve_mps(train_qubo, ns, 2; β=2., β_step=2, χ=2, threshold = 1e-14)
 
     #ground
     @test spins[1] == [1,-1,1,1,-1,1,1,1,1]
@@ -181,7 +189,9 @@ end
     end
     permuted_train_qubo = make_qubo1()
 
-    spins, objective = solve_mps(permuted_train_qubo, grid, 9, 16; β=2., β_step=2, χ=2, threshold = 1e-14)
+    ns = [Node_of_grid(i, permuted_train_qubo) for i in 1:get_system_size(permuted_train_qubo)]
+
+    spins, objective = solve_mps(permuted_train_qubo, ns, 16; β=2., β_step=2, χ=2, threshold = 1e-14)
 
     @test spins[1] == [1, -1, 1, 1, -1, 1, 1, 1, 1]
     @test objective[1] ≈ 0.30956452652382055
@@ -224,11 +234,13 @@ end
 
     println("number of β steps = ", β_step)
 
-    spins, _ = solve_mps(l_qubo, grid, 16, 10; β=β, β_step=β_step, χ=12, threshold = 1.e-8)
+    ns = [Node_of_grid(i, l_qubo) for i in 1:get_system_size(l_qubo)]
+
+    spins, _ = solve_mps(l_qubo, ns, 10; β=β, β_step=β_step, χ=12, threshold = 1.e-8)
 
     @test spins[1] == [1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1]
 
-    spins_exact, _ = solve_mps(l_qubo, grid, 16, 10; β=β, β_step=1, χ=12, threshold = 0.)
+    spins_exact, _ = solve_mps(l_qubo, ns, 10; β=β, β_step=1, χ=12, threshold = 0.)
 
     spins_peps, _ = solve(l_qubo, grid, 10; β = β, χ = 2, threshold = 1e-12)
 
@@ -239,4 +251,36 @@ end
         @test spins[k] == spins_peps[k]
     end
 end
+
+@testset "MPS vs PEPS larger QUBO" begin
+    function make_qubo_full()
+        qubo = [(1,1) 0.1; (1,2) 1.; (1,3) 1.; (1,4) 0.2; (2,2) -0.1; (2,3) 1.0; (2,4) 0.2]
+        qubo = vcat(qubo, [(3,3) 0.2; (3,4) 0.2; (4,4) -0.2])
+        [Qubo_el(qubo[i,1], qubo[i,2]) for i in 1:size(qubo, 1)]
+    end
+    l_qubo = make_qubo_full()
+
+    β = 0.5
+    β_step = 2
+
+    println("number of β steps = ", β_step)
+
+    ns = [Node_of_grid(i, l_qubo) for i in 1:get_system_size(l_qubo)]
+
+    spins, _ = solve_mps(l_qubo, ns, 4; β=β, β_step=β_step, χ=12, threshold = 1.e-8)
+
+    @test spins == [[1, 1, 1, 1], [-1, -1, -1, -1], [1, 1, 1, -1], [-1, -1, -1, 1]]
+
+    M = rand([-1.,-0.5,0.,0.5,1.], 64,64)
+    M = M*(M')
+
+    q = matrix2qubo_vec(M)
+
+    ns = [Node_of_grid(i, q) for i in 1:get_system_size(q)]
+
+    @time s, _ = solve_mps(q, ns, 4; β=β, β_step=β_step, χ=12, threshold = 1.e-8)
+    println(s[1])
+    println(s[2])
+    println(s[3])
+    println(s[4])
 end

src/tests/notation_tests.jl

@@ -2,6 +2,11 @@
 
     grid = [1 2 3; 4 5 6; 7 8 9]
 
+    b = Bond_with_other_node(1, [2,10], [3,11])
+    @test b.node == 1
+    @test b.spins1 == [2,10]
+    @test b.spins2 == [3,11]
+
     n = Node_of_grid(3, grid)
     @test n.i == 3
     @test n.spin_inds == [3]
@@ -10,6 +15,7 @@
     @test n.right == Array{Int64,1}[]
     @test n.up == Array{Int64,1}[]
     @test n.down == [[3,6]]
+    @test n.all_connections == [[3, 2], [3, 6]]
 
     n = Node_of_grid(5, grid)
     @test n.i == 5
@@ -61,6 +67,24 @@ function make_qubo0()
     [Qubo_el(qubo[i,1], qubo[i,2]) for i in 1:size(qubo, 1)]
 end
 
+@testset "axiliary on qubo" begin
+    qubo = make_qubo0()
+    n = Node_of_grid(1, qubo)
+    @test n.all_connections == [[1,2],[1,4]]
+
+    n = Node_of_grid(5, qubo)
+    @test n.all_connections == [[5, 2], [5, 6], [5, 4], [5, 8]]
+
+    n = Node_of_grid(9, qubo)
+    @test n.all_connections == [[9, 6], [9, 8]]
+
+    println(get_system_size(qubo) == 9)
+
+    M = rand(10,10)
+    matrix2qubo_vec(M)
+
+end
+
 @testset "test notation" begin
 
     @testset "Qubo_el type" begin

src/tests/test_full_graph.jl

@@ -0,0 +1,46 @@
+using NPZ
+using Plots
+using Test
+
+include("../notation.jl")
+include("../compression.jl")
+include("../peps_no_types.jl")
+include("../mps_implementation.jl")
+
+data = (x-> Array{Any}([parse(Int, x[1]), parse(Int, x[2]), parse(Float64, x[3])])).(split.(readlines(open("data/Qfile_tests.txt"))))
+
+function v2energy(M::Matrix{T}, v::Vector{Int}) where T <: AbstractFloat
+    d =  diag(M)
+    M = M .- diagm(d)
+
+    transpose(v)*M*v + transpose(v)*d
+end
+
+M = zeros(48,48)
+
+q_vec = Qubo_el{Float64}[]
+for d in data
+    i = ceil(Int, d[1])+1
+    j = ceil(Int, d[2])+1
+    if d[3] != 0.
+        push!(q_vec, Qubo_el((i,j), d[3]))
+        M[i,j] = d[3]
+    end
+end
+
+for i in 1:10
+    println(q_vec[i])
+end
+
+χ = 15
+β = 2.
+β_step = 2
+
+print("mps time  =  ")
+ns = [Node_of_grid(i, q_vec) for i in 1:get_system_size(q_vec)]
+@time spins_mps, objective_mps = solve_mps(q_vec, ns, 10; β=β, β_step=β_step, χ=χ, threshold = 1.e-8)
+
+println(spins_mps)
+for spin in spins_mps
+    println(v2energy(M, spin))
+end