Merge 3a14324 into 808459e

benchmarks/low_energy_approx_chimera.jl

@@ -11,7 +11,7 @@ using ArgParse
 using CSV
 using Test
 
-import SpinGlassPEPS: solve, solve_mps, M2graph, energy, binary2spins, ising_graph
+import SpinGlassPEPS: solve, M2graph, energy, binary2spins, ising_graph
 
 disable_logging(LogLevel(0))
 δH = 0.9

benchmarks/tests_grid_from_file.jl

@@ -10,7 +10,7 @@ using ArgParse
 using CSV
 using Test
 
-import SpinGlassPEPS: solve, solve_mps, M2graph
+import SpinGlassPEPS: solve, M2graph
 
 disable_logging(LogLevel(0))
 
@@ -193,7 +193,7 @@ for k in 1:examples
             @test states_given[i,:] == spins_spec[i]
         end
     end
-
+    #=
     ############ MPO - MPS #########
     χ = 15
 
@@ -220,4 +220,5 @@ for k in 1:examples
             @test states_given[i,:] == spins_mps[i]
         end
     end
+    =#
 end

benchmarks/tests_on_chimera.jl

@@ -10,7 +10,7 @@ using ArgParse
 using CSV
 using Test
 
-import SpinGlassPEPS: solve, solve_mps, M2graph, energy, binary2spins, ising_graph
+import SpinGlassPEPS: solve, M2graph, energy, binary2spins, ising_graph
 
 disable_logging(LogLevel(0))
 

src/SpinGlassPEPS.jl

@@ -9,7 +9,7 @@ module SpinGlassPEPS
     using Logging
     using StatsBase
     using Memoize, LRUCache
-    
+
     using DocStringExtensions
     const product = Iterators.product
 
@@ -27,7 +27,6 @@ module SpinGlassPEPS
     include("MPS_search.jl")
     include("notation.jl")
     include("peps_no_types.jl")
-    include("mps_implementation.jl")
 
     function __init__()
         @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" begin

src/contractions.jl

@@ -1,18 +1,18 @@
 export left_env, right_env, dot!
 
-# --------------------------- Conventions ------------------------ 
-#                                                                 
+# --------------------------- Conventions ------------------------
+#
 #      MPS          MPS*         MPO       left env      left env
 #       2            2            2           - 1          2 -
 #   1 - A - 3    1 - B - 3    1 - W - 3      L               R
 #                                 4           - 2          1 -
 # ---------------------------------------------------------------
 #
 
-function LinearAlgebra.dot(ψ::AbstractMPS, state::Union{Vector, NTuple}) 
+function LinearAlgebra.dot(ψ::AbstractMPS, state::Union{Vector, NTuple})
     C = I
-    
-    for (M, σ) ∈ zip(ψ, state)        
+
+    for (M, σ) ∈ zip(ψ, state)
         i = idx(σ)
         C = M[:, i, :]' * (C * M[:, i, :])
     end
@@ -26,12 +26,12 @@ function LinearAlgebra.dot(ϕ::AbstractMPS, ψ::AbstractMPS)
     for i ∈ eachindex(ψ)
         M = ψ[i]
         M̃ = conj(ϕ[i])
-        @tensor C[x, y] := M̃[β, σ, x] * C[β, α] * M[α, σ, y] order = (α, β, σ) 
+        @tensor C[x, y] := M̃[β, σ, x] * C[β, α] * M[α, σ, y] order = (α, β, σ)
     end
     C[]
 end
 
-function left_env(ϕ::AbstractMPS, ψ::AbstractMPS) 
+function left_env(ϕ::AbstractMPS, ψ::AbstractMPS)
     l = length(ψ)
     T = promote_type(eltype(ψ), eltype(ϕ))
 
@@ -54,7 +54,8 @@ end
     if l == 0
         L = [1.]
     else
-        m = idx(σ[l])
+        #m = idx(σ[l])
+        m = σ[l]
         L̃ = left_env(ϕ, σ[1:l-1])
         M = ϕ[l]
         @reduce L[x] := sum(α) L̃[α] * M[α, $m, x]
@@ -63,7 +64,7 @@ end
 end
 
 # NOT tested yet
-function right_env(ϕ::AbstractMPS, ψ::AbstractMPS) 
+function right_env(ϕ::AbstractMPS, ψ::AbstractMPS)
     L = length(ψ)
     T = promote_type(eltype(ψ), eltype(ϕ))
 
@@ -87,7 +88,8 @@ end
     if l == 0
         R = ones(1, 1)
     else
-        m = idx(σ[1])
+        #m = idx(σ[1])
+        m = σ[1]
         R̃ = right_env(ϕ, W, σ[2:l])
         M = ϕ[k-l+1]
         M̃ = W[k-l+1]
@@ -143,7 +145,7 @@ function LinearAlgebra.dot(O::AbstractMPO, ψ::AbstractMPS)
         W = O[i]
         M = ψ[i]
 
-        @reduce N[(x, a), σ, (y, b)] := sum(η) W[x, σ, y, η] * M[a, η, b]      
+        @reduce N[(x, a), σ, (y, b)] := sum(η) W[x, σ, y, η] * M[a, η, b]
         ϕ[i] = N
     end
     ϕ
@@ -155,7 +157,7 @@ function dot!(ψ::AbstractMPS, O::AbstractMPO)
         W = O[i]
         M = ψ[i]
 
-        @reduce N[(x, a), σ, (y, b)] := sum(η) W[x, σ, y, η] * M[a, η, b]      
+        @reduce N[(x, a), σ, (y, b)] := sum(η) W[x, σ, y, η] * M[a, η, b]
         ψ[i] = N
     end
 end
@@ -169,9 +171,9 @@ function LinearAlgebra.dot(O1::AbstractMPO, O2::AbstractMPO)
 
     for i ∈ 1:L
         W1 = O1[i]
-        W2 = O2[i]    
-        @reduce V[(x, a), σ, (y, b), η] := sum(γ) W1[x, σ, y, γ] * W2[a, γ, b, η]  
-                                        
+        W2 = O2[i]
+        @reduce V[(x, a), σ, (y, b), η] := sum(γ) W1[x, σ, y, γ] * W2[a, γ, b, η]
+
         O[i] = V
     end
     O

test/mps_tests.jl → src/detailed_mps/mps_tests.jl

@@ -1,7 +1,3 @@
-import SpinGlassPEPS: connections_for_mps, construct_mps
-import SpinGlassPEPS: contract4probability, compute_probs
-import SpinGlassPEPS: solve_mps, M2graph
-
 Random.seed!(1234)
 
 @testset "grouping of connections" begin

src/detailed_mps/runtests.jl

@@ -0,0 +1,18 @@
+using SpinGlassPEPS
+using LightGraphs
+using MetaGraphs
+using Random
+using Test
+using TensorOperations
+using LinearAlgebra
+using Statistics
+
+import SpinGlassPEPS: M2graph, Partial_sol, update_partial_solution, select_best_solutions
+
+include("../../test/test_helpers.jl")
+
+include("mps_implementation.jl")
+
+include("mps_tests.jl")
+include("tests_full_graph.jl")
+include("tests_on_data.jl")

src/peps_no_types.jl

@@ -29,103 +29,72 @@ function update_partial_solution(ps::Partial_sol{T}, s::Int, objective::T) where
     Partial_sol{T}(vcat(ps.spins, [s]), objective)
 end
 
-
-function project_spin_from_above(projector::Array{T,2}, spin::Int64, mps_el::Array{T,3}) where T <: Real
-    @reduce B[a, b] := sum(x) projector[$spin, x] * mps_el[a, x, b]
-    B
-end
-
-function project_spin_from_above(projector::Array{T,2}, spin::Int64, mpo_el::Array{T,4}) where T <: Real
-    @reduce B[a, d, c] := sum(x) projector[$spin, x] * mpo_el[a, x, c, d]
-    B
-end
-
-
 """
 ....
+    conditional_probabs(peps::PepsNetwork, ps::Partial_sol{T}, boundary_mps::MPS{T}, mpo::MPO{T}, peps_row::PEPSRow{T}) where T <: Number
 
 
-b_m - boundary mps
-p_r - peps row
+The contraction scheme
+
+mps - boundary mps
+T - 5 mode tensor
 
               physical .
-                         .   upper_right
-                           .  |      |
-    upper_left   from_left-- p_r3 --  p_r4 -- 1
-      |    |                 |      |
-1 --b_m1 --b_m2       --     b_m  -- b_m -- 1
+                         .       upper_right
+                           .  |              |
+    upper_left    from_left-- T --   ... -- mpo_peps.j_max -- 1
+      |    |                  |              |
+1 --mps_1-- mps_2 -- ...--   mps_i-- ... -- mps_peps.j_max -- 1
 """
-function conditional_probabs(peps::PepsNetwork, ps::Partial_sol{T}, boundary_mps::MPS{T}, peps_row::PEPSRow{T}) where T <: Number
-
+function conditional_probabs(peps::PepsNetwork, ps::Partial_sol{T}, boundary_mps::MPS{T}, mpo::MPO{T}, peps_row::PEPSRow{T}) where T <: Number
 
     j = length(ps.spins) + 1
     ng = peps.network_graph
     fg = ng.factor_graph
+    mpo1 = copy(mpo)
 
-
-    k = j % peps.j_max
-    if k == 0
-        k = peps.j_max
+    column = j % peps.j_max
+    if column == 0
+        column = peps.j_max
     end
     row = ceil(Int, j/peps.j_max)
 
-    mpo = MPO(peps_row)
-
     # set from above
     # not in last row
     if j > peps.j_max*(peps.i_max-1)
-        spin = [1 for _ in 1:k-1]
+        spin = [1 for _ in 1:column-1]
     else
-        spin = [ps.spins[i] for i in j-k+1:j-1]
-    end
-
-    proj_u = [projectors(fg, i, i+peps.j_max)[1] for i in j-k+1:j-1]
-
-
-    BB = [project_spin_from_above(proj_u[i], spin[i], boundary_mps[i]) for i in 1:k-1]
-
-    weight = ones(T, 1,1)
-    if k > 1
-        weight = prod(BB)
+        spin = [ps.spins[i] for i in j-column+1:j-1]
     end
+    proj_u = [projectors(fg, i, i+peps.j_max)[1] for i in j-column+1:j-1]
+    σ = Tuple(findall(l -> l == 1, proj_u[i][spin[i],:])[1] for i in 1:column-1)
+    w1 = left_env(boundary_mps, σ)
 
-    # set form left
-    if k == 1
-        spin = 1
-    else
-        spin = ps.spins[end]
-    end
-
-    proj_l, _, _ = projectors(fg, j-1, j)
-    @reduce A[d, a, b, c] := sum(x) proj_l[$spin, x] * peps_row[$k][x, a, b, c, d]
-
-    r = j-k+peps.j_max
+    # spinf for right env
+    r = j-column+peps.j_max
     if j <= peps.j_max
-        spin = [1 for _ in j:r]
+        spin_u = [1 for _ in j:r]
     else
-        spin = [ps.spins[i-peps.j_max] for i in j:r]
+        spin_u = [ps.spins[i-peps.j_max] for i in j:r]
     end
-
     proj_u = [projectors(fg, i-peps.j_max, i)[1] for i in j:r]
 
+    σ = Tuple(findall(l -> l == 1, proj_u[i][spin_u[i],:])[1] for i in 2:length(proj_u))
 
-    mpo = mpo[k+1:end]
-    mpo = MPO(vcat([A], mpo))
-
+    re1 = right_env(boundary_mps, mpo1, σ)
 
-    CC = [project_spin_from_above(proj_u[i], spin[i], mpo[i]) for i in 1:length(mpo)]
-
-    upper_mps = MPS(CC)
-
-    lower_mps = MPS(boundary_mps[k:end])
-
-    re = right_env(lower_mps, upper_mps)[1]
-
-    probs_unnormed = re*transpose(weight)
+    if column == 1
+        spin_l = 1
+    else
+        spin_l = ps.spins[end]
+    end
+    proj_l, _, _ = projectors(fg, j-1, j)
 
-    objective = probs_unnormed./sum(probs_unnormed)
+    ml = findall(l -> l == 1, proj_l[spin_l,:])[1]
+    mu = findall(u -> u == 1, proj_u[1][spin_u[1],:])[1]
+    @reduce P[σ] := sum(v, z, x, y) w1[v] * peps_row[$column][$ml, $mu, x, z, σ] * boundary_mps[$column][v, z, y]*re1[y,x]
 
-    dropdims(objective; dims = 2)
+    P./sum(P)
 end
 
 
@@ -201,6 +170,7 @@ function solve(peps::PepsNetwork, no_sols::Int = 2; β::T, χ::Int = typemax(Int
         @info "row of peps = " row
 
         peps_row = PEPSRow(peps, row)
+        mpo = MPO(peps_row)
 
         a = (row-1)*peps.j_max
 
@@ -214,7 +184,7 @@ function solve(peps::PepsNetwork, no_sols::Int = 2; β::T, χ::Int = typemax(Int
             partial_s = merge_dX(partial_s, dX, δH)
             for ps ∈ partial_s
 
-                objectives = conditional_probabs(peps, ps, boundary_mps[row], peps_row)
+                objectives = conditional_probabs(peps, ps, boundary_mps[row], mpo, peps_row)
 
                 for l ∈ eachindex(objectives)
                     new_objectives = ps.objective*objectives[l]