Low energy spectrum using PEPS (#73)

* rearange search jl * add new type AbstractGibbsNetwork * clean up search.jl * all test pass * start coding _δE * clean MPS_search * clean up a bit * add _δE * clean up * rm energy * correction in _δE * fix bug * clean up * add some comments * add correct update_energy * all test pass * rm idMPS, add improve MPS normalisation * improve energy * clean up * add test for low energy spectrum * fix bugs in search * add some corrections in search * MPS_search gives correct energy when max_states < N^2-1 * start dubegging prob * add tests * reavel problem with IdentityMPS * comment some test fot now * correct update_energy * Debug changes * clean up * start testing energy * try to find bug in energy * further attempts to improve energy * rzygi rules * Fixes * Add debugging code * add test for the gibbs state * clean up a bit * merge with origin/testing * adding missing bond-tensor in PEPS _contract * almost works! * the simplest example works * cleaning uo * add some tests to search_3.jl * clean up * some test stil do not pass * clean up * improve tests * fix the bug * clean up * the simples example works * improve tests * fix bugs in pathological instance * clean up searching examples * clean up * rm unnecassary function * clean up * pathological instance fo :NE works * clean up search * all tests pass, some broken * correct contraction example * setup not working example for today * Remove test deps * Finish merge * add some remarks * change order * search_2.jl works * ready to be merged * clean up * Delete create_sysimage( WTF was that. * incorporating minor improvments from @dexter * Add missing end * Update utils.jl Co-authored-by: Bartłomiej Gardas <[email protected]> Co-authored-by: annamariadziubyna <[email protected]> Co-authored-by: Konrad Jałowiecki <[email protected]> Co-authored-by: marekrams <[email protected]>
euro-hpc-pl · Mar 12, 2021 · 0822f62 · 0822f62
1 parent dc08660
commit 0822f62
Show file tree

Hide file tree

Showing 23 changed files with 438 additions and 240 deletions.
diff --git a/src/PEPS.jl b/src/PEPS.jl
@@ -8,15 +8,15 @@ const DEFAULT_CONTROL_PARAMS = Dict(
     "β" => 1.
 )
 
-struct PEPSNetwork
+struct PEPSNetwork <: AbstractGibbsNetwork
     size::NTuple{2, Int}
     map::Dict
     fg::MetaDiGraph
     nbrs::Dict
     origin::Symbol
     i_max::Int
     j_max::Int
-    β::Number
+    β::Number # TODO: get rid of this
     args::Dict{String, Number}
 
     function PEPSNetwork(
@@ -25,8 +25,8 @@ struct PEPSNetwork
         fg::MetaDiGraph,
         β::Number,
         origin::Symbol=:NW,
-        args_override::Dict{String, Number}=Dict{String, Number}()
-    )
+        args_override::Dict{String, T}=Dict{String, Number}()  # TODO: change String to Symbol
+    ) where T <: Number
         map, i_max, j_max = peps_indices(m, n, origin)
 
         # v => (l, u, r, d)
@@ -75,24 +75,25 @@ end
     else
         en = zeros(1, 1)
     end
-    exp.(-network.β .* en)
+    exp.(-network.β .* (en .- minimum(en)))
 end
 
-function PEPSRow(::Type{T}, peps::PEPSNetwork, i::Int) where {T <: Number}
-    ψ = PEPSRow(T, peps.j_max)
-
+function peps_tensor(::Type{T}, peps::PEPSNetwork, i::Int, j::Int) where {T <: Number}
     # generate tensors from projectors
-    for j ∈ 1:length(ψ)
-        ψ[j] = generate_tensor(peps, peps.map[i, j])
-    end
+    A = generate_tensor(peps, peps.map[i, j])
 
     # include energy
+    h = generate_tensor(peps, peps.map[i, j-1], peps.map[i, j])
+    v = generate_tensor(peps, peps.map[i-1, j], peps.map[i, j])
+    @tensor B[l, u, r, d, σ] := h[l, l̃] * v[u, ũ] * A[l̃, ũ, r, d, σ]
+    B
+end
+peps_tensor(peps::PEPSNetwork, i::Int, j::Int) = peps_tensor(Float64, peps, i, j)
+
+function PEPSRow(::Type{T}, peps::PEPSNetwork, i::Int) where {T <: Number}
+    ψ = PEPSRow(T, peps.j_max)
     for j ∈ 1:peps.j_max
-        A = ψ[j]
-        h = generate_tensor(peps, peps.map[i, j-1], peps.map[i, j])
-        v = generate_tensor(peps, peps.map[i-1, j], peps.map[i, j])
-        @tensor B[l, u, r, d, σ] := h[l, l̃] * v[u, ũ] * A[l̃, ũ, r, d, σ]
-        ψ[j] = B
+        ψ[j] = peps_tensor(T, peps, i, j)
     end
     ψ
 end
@@ -144,16 +145,13 @@ end
 function contract_network(
     peps::PEPSNetwork,
     config::Dict{Int, Int} = Dict{Int, Int}(),
-    )
+)
     ψ = MPS(peps, 1, config)
-    prod(dropdims(ψ))[]
+    prod(dropindices(ψ))[]
 end
 
-@inline function _get_coordinates(
-    peps::PEPSNetwork,
-    k::Int
-    )
-    ceil(k / peps.j_max), (k - 1) % peps.j_max + 1
+@inline function get_coordinates(peps::PEPSNetwork, k::Int)
+    ceil(Int, k / peps.j_max), (k - 1) % peps.j_max + 1
 end
 
 function generate_boundary(fg::MetaDiGraph, v::Int, w::Int, state::Int)
@@ -163,7 +161,8 @@ function generate_boundary(fg::MetaDiGraph, v::Int, w::Int, state::Int)
     findfirst(x -> x > 0, pv[state, :])
 end
 
-function generate_boundary(peps::PEPSNetwork, v::Vector{Int}, i::Int, j::Int)
+function generate_boundary(peps::PEPSNetwork, v::Vector{Int}, w::NTuple{2, Int})
+    i, j = w
     ∂v = zeros(Int, peps.j_max + 1)
 
     # on the left below
@@ -172,45 +171,33 @@ function generate_boundary(peps::PEPSNetwork, v::Vector{Int}, i::Int, j::Int)
             peps.fg,
             peps.map[i, k],
             peps.map[i+1, k],
-            _get_local_state(peps, v, i, k))
+            _get_local_state(peps, v, (i, k))
+        )
     end
 
     # on the left at the current row
     ∂v[j] = generate_boundary(
         peps.fg,
         peps.map[i, j-1],
         peps.map[i, j],
-        _get_local_state(peps, v, i, j-1))
+        _get_local_state(peps, v, (i, j-1))
+    )
 
     # on the right above
     for k ∈ j:peps.j_max
         ∂v[k+1] = generate_boundary(
             peps.fg,
             peps.map[i-1, k],
             peps.map[i, k],
-            _get_local_state(peps, v, i-1, k))
+            _get_local_state(peps, v, (i-1, k))
+        )
     end
     ∂v
 end
 
-function _get_local_state(peps::PEPSNetwork, v::Vector{Int}, i::Int, j::Int)
-    k = j + peps.j_max * (i - 1)
-    0 < k <= lenght(v) ? v[k] : 1
-end
-
-# function generate_boundary(peps::PEPSNetwork, v::Vector{Int})
-#     i, j = _get_coordinates(peps, length(v)+1)
-#     generate_boundary(peps, v, i, j)
-# end
-
-function _contract(
-    A::Array{T, 5}, M::Array{T, 3}, L::Vector{T}, R::Matrix{T}, ∂v::Vector{Int}
-) where {T <: Number}
-    l, u = ∂v
-    @cast Ã[r, d, σ] := A[$l, $u, r, d, σ]
-    @tensor prob[σ] := L[x] * M[x, d, y] *
-                       Ã[r, d, σ] * R[y, r] order = (x, d, r, y)
-    prob
+function _get_local_state(peps::PEPSNetwork, σ::Vector{Int}, w::NTuple{2, Int})
+    k = w[2] + peps.j_max * (w[1] - 1)
+    0 < k <= length(σ) ? σ[k] : 1
 end
 
 function _normalize_probability(prob::Vector{T}) where {T <: Number}
@@ -223,24 +210,52 @@ function conditional_probability(
     peps::PEPSNetwork,
     v::Vector{Int},
     )
-    i, j = _get_coordinates(peps, length(v)+1)
-    ∂v = generate_boundary(peps, v, i, j)
+    i, j = get_coordinates(peps, length(v)+1)
+    ∂v = generate_boundary(peps, v, (i, j))
 
     W = MPO(peps, i)
     ψ = MPS(peps, i+1)
 
     L = left_env(ψ, ∂v[1:j-1])
     R = right_env(ψ, W, ∂v[j+2:peps.j_max+1])
-    A = generate_tensor(peps, i, j)
+    A = peps_tensor(peps, i, j)
 
     l, u = ∂v[j:j+1]
     M = ψ[j]
-    Ã[:, :, :] = A[l, u, :, :, :]
+    Ã = A[l, u, :, :, :]
     @tensor prob[σ] := L[x] * M[x, d, y] *
                        Ã[r, d, σ] * R[y, r] order = (x, d, r, y)
     _normalize_probability(prob)
 end
 
+function bond_energy(fg::MetaDiGraph, u::Int, v::Int, σ::Int)
+    if has_edge(fg, u, v)
+        pu, en, pv = get_prop.(Ref(fg), u, v, (:pl, :en, :pr))
+        energies = (pu * (en * pv[:, σ:σ]))'
+    elseif has_edge(fg, v, u)
+        pv, en, pu = get_prop.(Ref(fg), v, u, (:pl, :en, :pr))
+        energies = (pv[σ:σ, :] * en) * pu
+    else
+        energies = zeros(get_prop(fg, u, :loc_dim))
+    end
+    vec(energies)
+end
+
+function update_energy(
+    network::AbstractGibbsNetwork,
+    σ::Vector{Int},
+    )
+    i, j = get_coordinates(network, length(σ)+1)
+
+    σkj = _get_local_state(network, σ, (i-1, j))
+    σil = _get_local_state(network, σ, (i, j-1))
+
+    bond_energy(network.fg, network.map[i, j], network.map[i, j-1], σil) +
+    bond_energy(network.fg, network.map[i, j], network.map[i-1, j], σkj) +
+    get_prop(network.fg, network.map[i, j], :loc_en)
+end
+
+#TODO: translate this into rotations and reflections
 function peps_indices(m::Int, n::Int, origin::Symbol=:NW)
     @assert origin ∈ (:NW, :WN, :NE, :EN, :SE, :ES, :SW, :WS)
 

diff --git a/src/SpinGlassPEPS.jl b/src/SpinGlassPEPS.jl
@@ -22,6 +22,7 @@ module SpinGlassPEPS
     include("ising.jl")
     include("exact.jl")
     include("factor.jl")
+    include("search.jl")
     include("PEPS.jl")
     include("MPS_search.jl")
 

diff --git a/src/base.jl b/src/base.jl
@@ -1,6 +1,6 @@
 export bond_dimension, is_left_normalized, is_right_normalized
 export verify_bonds, verify_physical_dims, tensor, rank, physical_dim
-export State, idMPS
+export State, dropindices
 
 const State = Union{Vector, NTuple}
 
@@ -50,9 +50,7 @@ end
 @inline MPS(A::AbstractArray, ::Val{:right}, Dcut::Int, args...) = _left_sweep_SVD(MPS, A, Dcut, args...)
 @inline MPS(A::AbstractArray, ::Val{:left}, Dcut::Int, args...) = _right_sweep_SVD(MPS, A, Dcut, args...)
 
-@inline Base.dropdims(ψ::MPS, i::Int) = (dropdims(A, dims=i) for A ∈ ψ)
-@inline Base.dropdims(ψ::MPS) = Base.dropdims(ψ, 2)
-
+@inline dropindices(ψ::AbstractMPS, i::Int=2) = (dropdims(A, dims=i) for A ∈ ψ)
 
 function MPS(states::Vector{Vector{T}}) where {T <: Number}
     state_arrays = [reshape(copy(v), (1, length(v), 1)) for v ∈ states]

diff --git a/src/compressions.jl b/src/compressions.jl
@@ -50,7 +50,7 @@ function _right_sweep_SVD!(ψ::AbstractMPS, Dcut::Int=typemax(Int))
 
     for i ∈ eachindex(ψ)
         A = ψ[i]
-        C = Diagonal(Σ) * V'
+        C = (Diagonal(Σ) ./ Σ[1]) * V'
 
         # attach
         @tensor M[x, σ, y] := C[x, α] * A[α, σ, y]
@@ -64,14 +64,15 @@ function _right_sweep_SVD!(ψ::AbstractMPS, Dcut::Int=typemax(Int))
         @cast A[x, σ, y] |= U[(x, σ), y] (σ:d)
         ψ[i] = A
     end
+    ψ[end] *= tr(V)
 end
 
 function _left_sweep_SVD!(ψ::AbstractMPS, Dcut::Int=typemax(Int))
     Σ = U = ones(eltype(ψ), 1, 1)
 
     for i ∈ length(ψ):-1:1
         B = ψ[i]
-        C = U * Diagonal(Σ)
+        C = U * (Diagonal(Σ) ./ Σ[1])
 
         # attach
         @tensor M[x, σ, y]   := B[x, σ, α] * C[α, y]
@@ -85,6 +86,7 @@ function _left_sweep_SVD!(ψ::AbstractMPS, Dcut::Int=typemax(Int))
         @cast B[x, σ, y] |= V'[x, (σ, y)] (σ:d)
         ψ[i] = B
     end
+    ψ[1] *= tr(U)
 end
 
 function _left_sweep_var!!(ϕ::AbstractMPS, env::Vector{<:AbstractMatrix}, ψ::AbstractMPS, Dcut::Int)

diff --git a/src/contractions.jl b/src/contractions.jl
@@ -1,7 +1,7 @@
 export left_env, right_env, dot!
 
 # --------------------------- Conventions -----------------------
-#                                                                 
+#
 #      MPS          MPS*         MPO       left env     right env
 #       2            2            2           - 1          2 -
 #   1 - A - 3    1 - B - 3    1 - W - 3      L               R
@@ -85,7 +85,8 @@ end
 
 @memoize function right_env(ϕ::AbstractMPS{T}, W::AbstractMPO{T}, σ::Union{Vector, NTuple}) where {T}
     l = length(σ)
-    k = length(ϕ)
+    #k = length(ϕ)
+    k = length(W)
     if l == 0
         R = similar(ϕ[1], T, (1, 1))
         R[1, 1] = one(T)

diff --git a/src/exact.jl b/src/exact.jl
@@ -43,7 +43,7 @@ function brute_force(ig::MetaGraph; sorted=true, num_states::Int=1)
         perm = partialsortperm(vec(energies), 1:num_states)
         return Spectrum(energies[perm], σ[perm])
     else
-        return Spectrum(energies, σ)
+        return Spectrum(energies[1:num_states], σ[1:num_states])
     end
 end
 

diff --git a/src/factor.jl b/src/factor.jl
@@ -58,6 +58,7 @@ function factor_graph(
         sp = spectrum(cl, num_states=get(num_states_cl, v, basis_size(cl)))
         set_prop!(fg, v, :spectrum, sp)
         set_prop!(fg, v, :loc_en, vec(sp.energies))
+        set_prop!(fg, v, :loc_dim, length(vec(sp.energies)))
     end
 
     for i ∈ 1:L, j ∈ i+1:L
@@ -82,19 +83,6 @@ function factor_graph(
     fg
 end
 
-
-# TODO: eradicate it
-function projectors(fg::MetaDiGraph, i::Int, j::Int)
-    if has_edge(fg, i, j)
-        p1, en, p2 = get_prop(fg, i, j, :split)
-    elseif has_edge(fg, j, i)
-        p2, en, p1 = get_prop(fg, j, i, :split)
-    else
-        p1 = en = p2 = ones(1, 1)
-    end
-    p1, en, p2
-end
-
 function rank_reveal(energy, order=:PE)
     @assert order ∈ (:PE, :EP)
     dim = order == :PE ? 1 : 2

diff --git a/src/ising.jl b/src/ising.jl
@@ -5,8 +5,8 @@ export Spectrum, cluster, rank, nodes, basis_size
 const Instance = Union{String, Dict}
 
 struct Spectrum
-    energies::Array{<:Number}
-    states::Array{Vector{<:Number}}
+    energies::Vector{Float64}
+    states::Vector{Vector{Int}}
 end
 
 unique_nodes(ising_tuples) = sort(collect(Set(Iterators.flatten((i, j) for (i, j, _) ∈ ising_tuples))))

diff --git a/src/lattice.jl b/src/lattice.jl
@@ -1,6 +1,6 @@
-export chimera_to_square_lattice
+export super_square_lattice
 
- function chimera_to_square_lattice(size::NTuple{5, Int})
+ function super_square_lattice(size::NTuple{5, Int})
     m, um, n, un, t = size
     new = LinearIndices((1:n, 1:m))
     old = LinearIndices((1:t, 1:un, 1:n, 1:um, 1:m))
@@ -11,7 +11,7 @@ export chimera_to_square_lattice
     )
 end
 
-function chimera_to_square_lattice(size::NTuple{3, Int})
+function super_square_lattice(size::NTuple{3, Int})
     m, n, t = size
-    chimera_to_square_lattice((m, 1, n, 1, t))
+    super_square_lattice((m, 1, n, 1, t))
 end