From 4a889a89c47b2608a9c2554212391ca71d2f4c7c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Konrad=20Ja=C5=82owiecki?= Date: Tue, 9 Mar 2021 23:33:22 +0100 Subject: [PATCH] Machete kills network graph --- src/PEPS.jl | 146 +++++++++++++++++++++++-------------------- src/SpinGlassPEPS.jl | 3 - src/factor.jl | 2 +- src/ising.jl | 8 +-- src/network.jl | 55 ---------------- test/factor.jl | 5 +- 6 files changed, 83 insertions(+), 136 deletions(-) delete mode 100644 src/network.jl diff --git a/src/PEPS.jl b/src/PEPS.jl index f45c6e8..ea4004b 100644 --- a/src/PEPS.jl +++ b/src/PEPS.jl @@ -8,7 +8,7 @@ const DEFAULT_CONTROL_PARAMS = Dict( "β" => 1. ) -mutable struct PEPSNetwork +struct PEPSNetwork size::NTuple{2, Int} map::Dict fg::MetaDiGraph @@ -26,54 +26,71 @@ mutable struct PEPSNetwork β::Number, origin::Symbol=:NW, args_override::Dict{String, Number}=Dict{String, Number}() + ) + map, i_max, j_max = peps_indices(m, n, origin) + + # v => (l, u, r, d) + nbrs = Dict( + map[i, j] => (map[i, j-1], map[i-1, j], map[i, j+1], map[i+1, j]) + for i ∈ 1:i_max, j ∈ 1:j_max ) - pn = new((m, n)) - pn.map, pn.i_max, pn.j_max = peps_indices(m, n, origin) + args = merge(DEFAULT_CONTROL_PARAMS, args_override) + pn = new((m, n), map, fg, nbrs, origin, i_max, j_max, β, args) + end +end - nbrs = Dict() - for i ∈ 1:pn.i_max, j ∈ 1:pn.j_max - # v => (l, u, r, d) - 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.fg = fg - pn.nbrs = nbrs - pn.β = β - pn.args = merge(DEFAULT_CONTROL_PARAMS, args_override) - pn +function _get_projector(fg::MetaDiGraph, v::Int, w::Int) + if has_edge(fg, w, v) + get_prop(fg, w, v, :pr)' + elseif has_edge(fg, v, w) + get_prop(fg, v, w, :pl) + else + loc_dim = length(get_prop(fg, v, :loc_en)) + ones(loc_dim, 1) end end -generate_tensor(pn::PEPSNetwork, - m::NTuple{2,Int}, - ) = generate_tensor(pn, pn.map[m]) +@memoize function generate_tensor(network::PEPSNetwork, v::Int) + # TODO: does this require full network, or can we pass only fg? + loc_exp = exp.(-network.β .* get_prop(network.fg, v, :loc_en)) -generate_tensor(pn::PEPSNetwork, - m::NTuple{2, Int}, - n::NTuple{2, Int}, - ) = generate_tensor(pn, pn.map[m], pn.map[n]) + dim = zeros(Int, length(network.nbrs[v])) + @cast A[_, i] := loc_exp[i] -generate_boundary(pn::PEPSNetwork, - m::NTuple{2, Int}, - n::NTuple{2, Int}, - σ::Int, - ) = generate_boundary(pn, pn.map[m], pn.map[n], σ) + for (j, w) ∈ enumerate(network.nbrs[v]) + pv = _get_projector(network.fg, v, w) + @cast A[(c, γ), σ] |= A[c, σ] * pv[σ, γ] + dim[j] = size(pv, 2) + end + reshape(A, dim..., :) +end + +@memoize function generate_tensor(network::PEPSNetwork, v::Int, w::Int) + fg = network.fg + if has_edge(fg, w, v) + en = get_prop(fg, w, v, :en)' + elseif has_edge(fg, v, w) + en = get_prop(fg, v, w, :en) + else + en = zeros(1, 1) + end + exp.(-network.β .* en) +end function PEPSRow(::Type{T}, peps::PEPSNetwork, i::Int) where {T <: Number} ψ = PEPSRow(T, peps.j_max) # generate tensors from projectors for j ∈ 1:length(ψ) - ψ[j] = generate_tensor(peps, (i, j)) + ψ[j] = generate_tensor(peps, peps.map[i, j]) end # include energy for j ∈ 1:peps.j_max A = ψ[j] - h = generate_tensor(peps, (i, j-1), (i, j)) - v = generate_tensor(peps, (i-1, j), (i, 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, ũ] * A[l̃, ũ, r, d, σ] ψ[j] = B end @@ -101,6 +118,7 @@ function MPO(::Type{T}, end W end + MPO(peps::PEPSNetwork, i::Int, config::Dict{Int, Int} = Dict{Int, Int}() @@ -138,68 +156,56 @@ end ceil(k / peps.j_max), (k - 1) % peps.j_max + 1 end -@inline function _get_local_state( - peps::PEPSNetwork, - v::Vector{Int}, - i::Int, - j::Int, - ) - k = j + peps.j_max * (i - 1) - if k > length(v) || k <= 0 return 1 end - v[k] +function generate_boundary(fg::MetaDiGraph, v::Int, w::Int, state::Int) + if v ∉ vertices(fg) return 1 end + loc_dim = length(get_prop(fg, v, :loc_en)) + pv = _get_projector(fg, v, w) + 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}, i::Int, j::Int) ∂v = zeros(Int, peps.j_max + 1) # on the left below for k ∈ 1:j-1 ∂v[k] = generate_boundary( - peps.network_graph, - (i, k), - (i+1, k), + peps.fg, + peps.map[i, k], + peps.map[i+1, k], _get_local_state(peps, v, i, k)) end # on the left at the current row ∂v[j] = generate_boundary( - peps.network_graph, - (i, j-1), - (i, j), + peps.fg, + peps.map[i, j-1], + peps.map[i, j], _get_local_state(peps, v, i, j-1)) # on the right above for k ∈ j:peps.j_max ∂v[k+1] = generate_boundary( - peps.network_graph, - (i-1, k), - (i, k), + peps.fg, + peps.map[i-1, k], + peps.map[i, k], _get_local_state(peps, v, i-1, k)) end ∂v end -function generate_boundary( - peps::PEPSNetwork, - v::Vector{Int}, - ) - i, j = _get_coordinates(peps, length(v)+1) - generate_boundary(peps, v, i, j) +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 -@inline function _contract( - A::Array{T, 5}, - M::Array{T, 3}, - L::Vector{T}, - R::Matrix{T}, - ∂v::Vector{Int}, - ) where {T <: Number} +# 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 Ã[r, d, σ] := A[$l, $u, r, d, σ] @tensor prob[σ] := L[x] * M[x, d, y] * @@ -227,7 +233,11 @@ function conditional_probability( R = right_env(ψ, W, ∂v[j+2:peps.j_max+1]) A = generate_tensor(peps, i, j) - prob = _contract(A, ψ[j], L, R, ∂v[j:j+1]) + l, u = ∂v[j:j+1] + M = ψ[j] + Ã[:, :, :] = A[l, u, :, :, :] + @tensor prob[σ] := L[x] * M[x, d, y] * + Ã[r, d, σ] * R[y, r] order = (x, d, r, y) _normalize_probability(prob) end diff --git a/src/SpinGlassPEPS.jl b/src/SpinGlassPEPS.jl index 49dfa26..e456582 100644 --- a/src/SpinGlassPEPS.jl +++ b/src/SpinGlassPEPS.jl @@ -23,10 +23,7 @@ module SpinGlassPEPS include("exact.jl") include("factor.jl") include("PEPS.jl") - include("network.jl") include("MPS_search.jl") - include("notation.jl") - include("peps_no_types.jl") function __init__() @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" begin diff --git a/src/factor.jl b/src/factor.jl index 8245779..511c2d1 100644 --- a/src/factor.jl +++ b/src/factor.jl @@ -32,7 +32,7 @@ function factor_graph( energy::Function=energy, spectrum::Function=full_spectrum, cluster_assignment_rule::Dict{Int, Int} # e.g. square lattice - ) +) ns = Dict(i => num_states_cl for i ∈ Set(values(cluster_assignment_rule))) factor_graph( ig, diff --git a/src/ising.jl b/src/ising.jl index 46dffb0..a3e390c 100644 --- a/src/ising.jl +++ b/src/ising.jl @@ -3,8 +3,6 @@ export energy, rank_vec export Spectrum, cluster, rank, nodes, basis_size const Instance = Union{String, Dict} -const SimpleEdge = LightGraphs.SimpleGraphs.SimpleEdge -const EdgeIter = Union{LightGraphs.SimpleGraphs.SimpleEdgeIter, Base.Iterators.Filter, Array} struct Spectrum energies::Array{<:Number} @@ -85,11 +83,7 @@ function cluster(ig::MetaGraph, verts) rank = getindex.(Ref(get_prop(ig, :rank)), vmap) J = get_prop(ig, :J)[vmap, vmap] - set_prop!(sub_ig, :rank, rank) - set_prop!(sub_ig, :J, J) - set_prop!(sub_ig, :h, h) - set_prop!(sub_ig, :vmap, vmap) - + set_props!(sub_ig, Dict(:rank => rank, :J => J, :h => h, :vmap => vmap)) sub_ig end diff --git a/src/network.jl b/src/network.jl deleted file mode 100644 index 12ab4bc..0000000 --- a/src/network.jl +++ /dev/null @@ -1,55 +0,0 @@ -export generate_boundary, generate_tensor - - -function _get_projector( - fg::MetaDiGraph, - v::Int, - w::Int - ) - if has_edge(fg, w, v) - get_prop(fg, w, v, :pr)' - elseif has_edge(fg, v, w) - get_prop(fg, v, w, :pl) - else - loc_dim = length(get_prop(fg, v, :loc_en)) - ones(loc_dim, 1) - end -end - -@memoize function generate_tensor(network::PEPSNetwork, v::Int) - loc_exp = exp.(-network.β .* get_prop(network.fg, v, :loc_en)) - - dim = zeros(Int, length(network.nbrs[v])) - @cast A[_, i] := loc_exp[i] - - for (j, w) ∈ enumerate(network.nbrs[v]) - pv = _get_projector(network.fg, v, w) - @cast A[(c, γ), σ] |= A[c, σ] * pv[σ, γ] - dim[j] = size(pv, 2) - end - reshape(A, dim..., :) -end - -@memoize function generate_tensor(network::PEPSNetwork, v::Int, w::Int) - fg = network.fg - if has_edge(fg, w, v) - en = get_prop(fg, w, v, :en)' - elseif has_edge(fg, v, w) - en = get_prop(fg, v, w, :en) - else - en = zeros(1, 1) - end - exp.(-network.β .* en) -end - -function generate_boundary( - fg::MetaDiGraph, - v::Int, - w::Int, - state::Int - ) - if v ∉ vertices(fg) return 1 end - loc_dim = length(get_prop(fg, v, :loc_en)) - pv = _get_projector(fg, v, w) - findfirst(x -> x > 0, pv[state, :]) -end diff --git a/test/factor.jl b/test/factor.jl index 20ed63d..b524810 100644 --- a/test/factor.jl +++ b/test/factor.jl @@ -166,8 +166,9 @@ end @testset "each edge comprises expected bunch of edges from source Ising graph" begin for e ∈ edges(fg) outer_edges = get_prop(fg, e, :outer_edges) - println(collect(outer_edges)) - println(cedges[Tuple(e)]) + # println(collect(outer_edges)) + # println(cedges[Tuple(e)]) + # Note: this test is ok if we translate edges correctly. # TODO: fix this by translating from nodes to graph coordinates # @test issetequal(cedges[Tuple(e)], collect(outer_edges)) end