From eec0e2cc786ece00561b10cdd1b79e932ebe4b9c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C5=81ukasz=20Pawela?= <3093117+lpawela@users.noreply.github.com> Date: Wed, 14 Apr 2021 12:26:31 +0200 Subject: [PATCH] Revert "[WIP] new network interface" (#5) * Revert "[WIP] new network interface (#4)" This reverts commit 0d891e13a4e05d228e29f6a4b3521ac260e114ce. * Update Project.toml --- Project.toml | 1 - src/MPS_search.jl | 71 +++++++-- src/PEPS.jl | 297 +++++++++++++++++-------------------- src/SpinGlassEngine.jl | 2 - src/attic.jl | 54 ------- src/network_interface.jl | 97 ------------ src/network_operations.jl | 84 ----------- src/search.jl | 26 ++-- test/PEPS.jl | 13 +- test/contract.jl | 21 +-- test/network_operations.jl | 126 ---------------- test/runtests.jl | 3 +- test/search.jl | 46 +++--- test/search_2.jl | 18 +-- 14 files changed, 255 insertions(+), 604 deletions(-) delete mode 100644 src/attic.jl delete mode 100644 src/network_interface.jl delete mode 100644 src/network_operations.jl delete mode 100644 test/network_operations.jl diff --git a/Project.toml b/Project.toml index b8907d1a..674d5141 100644 --- a/Project.toml +++ b/Project.toml @@ -4,7 +4,6 @@ version = "0.1.0" [deps] DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -LabelledGraphs = "605abd48-4d17-4660-b914-d4df33194460" LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Memoize = "c03570c3-d221-55d1-a50c-7939bbd78826" diff --git a/src/MPS_search.jl b/src/MPS_search.jl index 9099a17b..a8eea2fd 100644 --- a/src/MPS_search.jl +++ b/src/MPS_search.jl @@ -9,10 +9,63 @@ struct MPSControl dβ::Number end +_make_left_env(ψ::AbstractMPS, k::Int) = ones(eltype(ψ), 1, 1, k) +_make_left_env_new(ψ::AbstractMPS, k::Int) = ones(eltype(ψ), 1, k) +_make_LL(ψ::AbstractMPS, b::Int, k::Int, d::Int) = zeros(eltype(ψ), b, b, k, d) +_make_LL_new(ψ::AbstractMPS, b::Int, k::Int, d::Int) = zeros(eltype(ψ), b, k, d) -_make_left_env(ψ::AbstractMPS, k::Int) = ones(eltype(ψ), 1, k) +# ρ needs to be ∈ the right canonical form +# function solve(ψ::AbstractMPS, keep::Int) +# @assert keep > 0 "Number of states has to be > 0" +# T = eltype(ψ) + +# keep_extra = keep +# pCut = prob = 0. +# k = 1 + +# if keep < prod(rank(ψ)) +# keep_extra += 1 +# end + +# states = fill([], 1, k) +# left_env = _make_left_env(ψ, k) + +# for (i, M) ∈ enumerate(ψ) +# _, d, b = size(M) + +# pdo = zeros(T, k, d) +# LL = _make_LL(ψ, b, k, d) +# config = zeros(Int, i, k, d) + +# for j ∈ 1:k +# L = left_env[:, :, j] + +# for σ ∈ local_basis(d) +# m = idx(σ) +# LL[:, :, j, m] = M[:, m, :]' * (L * M[:, m, :]) +# pdo[j, m] = tr(LL[:, :, j, m]) +# config[:, j, m] = vcat(states[:, j]..., σ) +# end +# end + +# perm = collect(1: k * d) +# k = min(k * d, keep_extra) + +# if k >= keep_extra +# partialsortperm!(perm, vec(pdo), 1:k, rev=true) +# prob = vec(pdo)[perm] +# pCut < last(prob) ? pCut = last(prob) : () +# end + +# @cast A[α, β, (l, d)] |= LL[α, β, l, d] +# left_env = A[:, :, perm] + +# @cast B[α, (l, d)] |= config[α, l, d] +# states = B[:, perm] +# end +# states[:, 1:keep], prob[1:keep], pCut +# end -_make_LL(ψ::AbstractMPS, b::Int, k::Int, d::Int) = zeros(eltype(ψ), b, k, d) # ψ needs to be ∈ the right canonical form function solve(ψ::AbstractMPS, keep::Int) @@ -30,14 +83,14 @@ function solve(ψ::AbstractMPS, keep::Int) lprob = zeros(T, k) states = fill([], 1, k) - left_env = _make_left_env(ψ, k) + left_env = _make_left_env_new(ψ, k) for (i, M) ∈ enumerate(ψ) _, d, b = size(M) pdo = ones(T, k, d) lpdo = zeros(T, k, d) - LL = _make_LL(ψ, b, k, d) + LL = _make_LL_new(ψ, b, k, d) config = zeros(Int, i, k, d) for j ∈ 1:k @@ -73,7 +126,7 @@ function solve(ψ::AbstractMPS, keep::Int) states', lprob, lpCut end -function _apply_bias!(ψ::AbstractMPS, ig::LabelledGraph, dβ::Number, i::Int) +function _apply_bias!(ψ::AbstractMPS, ig::MetaGraph, dβ::Number, i::Int) M = ψ[i] d = size(M, 2) @@ -84,7 +137,7 @@ function _apply_bias!(ψ::AbstractMPS, ig::LabelledGraph, dβ::Number, i::Int) ψ[i] = M end -function _apply_exponent!(ψ::AbstractMPS, ig::LabelledGraph, dβ::Number, i::Int, j::Int, last::Int) +function _apply_exponent!(ψ::AbstractMPS, ig::MetaGraph, dβ::Number, i::Int, j::Int, last::Int) M = ψ[j] D = typeof(M).name.wrapper(I(physical_dim(ψ, i))) @@ -133,7 +186,7 @@ end _holes(l::Int, nbrs::Vector) = setdiff(l+1 : last(nbrs), nbrs) -function _apply_layer_of_gates(ig::LabelledGraph, ρ::AbstractMPS, control::MPSControl, dβ::Number) +function _apply_layer_of_gates(ig::MetaGraph, ρ::AbstractMPS, control::MPSControl, dβ::Number) L = nv(ig) Dcut = control.max_bond tol = control.var_ϵ @@ -168,7 +221,7 @@ function _apply_layer_of_gates(ig::LabelledGraph, ρ::AbstractMPS, control::MPSC ρ end -function SpinGlassTensors.MPS(ig::LabelledGraph, control::MPSControl) +function SpinGlassTensors.MPS(ig::MetaGraph, control::MPSControl) Dcut = control.max_bond tol = control.var_ϵ @@ -187,7 +240,7 @@ function SpinGlassTensors.MPS(ig::LabelledGraph, control::MPSControl) ρ end -function SpinGlassTensors.MPS(ig::LabelledGraph, control::MPSControl, type::Symbol) +function SpinGlassTensors.MPS(ig::MetaGraph, control::MPSControl, type::Symbol) L = nv(ig) Dcut = control.max_bond tol = control.var_ϵ diff --git a/src/PEPS.jl b/src/PEPS.jl index da326fae..dfea81e7 100644 --- a/src/PEPS.jl +++ b/src/PEPS.jl @@ -8,138 +8,108 @@ const DEFAULT_CONTROL_PARAMS = Dict( "β" => 1. ) -# struct PEPSNetwork <: AbstractGibbsNetwork -# size::NTuple{2, Int} -# map::Dict -# fg::MetaDiGraph -# nbrs::Dict -# origin::Symbol -# i_max::Int -# j_max::Int -# β::Number # TODO: get rid of this -# args::Dict{String, Number} - -# function PEPSNetwork( -# m::Int, -# n::Int, -# fg::MetaDiGraph, -# β::Number, -# origin::Symbol=:NW, -# 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) -# 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 -# ) - -# args = merge(DEFAULT_CONTROL_PARAMS, args_override) -# pn = new((m, n), map, fg, nbrs, origin, i_max, j_max, β, args) -# end -# end - -function peps_lattice(m, n) - labels = [(i, j) for j ∈ 1:n for i ∈ 1:m] - LabelledGraph(labels, grid((m, n))) +struct PEPSNetwork <: AbstractGibbsNetwork + size::NTuple{2, Int} + map::Dict + fg::MetaDiGraph + nbrs::Dict + origin::Symbol + i_max::Int + j_max::Int + β::Number # TODO: get rid of this + args::Dict{String, Number} + + function PEPSNetwork( + m::Int, + n::Int, + fg::MetaDiGraph, + β::Number, + origin::Symbol=:NW, + 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) + 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 + ) + + args = merge(DEFAULT_CONTROL_PARAMS, args_override) + pn = new((m, n), map, fg, nbrs, origin, i_max, j_max, β, args) + end end - -struct PEPSNetwork <: AbstractGibbsNetwork{NTuple{2, Int}, NTuple{2, Int}} - factor_graph::LabelledGraph{T, NTuple{2, Int}} where T - network_graph::LabelledGraph{S, NTuple{2, Int}} where S - vertex_map::Function - m::Int - n::Int - nrows::Int - ncols::Int - - function PEPSNetwork(m::Int, n::Int, factor_graph, transformation::LatticeTransformation) - vmap = vertex_map(transformation, m, n) - ng = peps_lattice(m, n) - nrows, ncols = transformation.flips_dimensions ? (n, m) : (m, n) - if !is_compatible(factor_graph, ng, vmap) - throw(ArgumentError("Factor graph not compatible with given network.")) - end - new(factor_graph, ng, vmap, m, n, nrows, ncols) +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) + # TODO: does this require full network, or can we pass only fg? + loc_exp = exp.(-network.β .* get_prop(network.fg, v, :loc_en)) + dim = zeros(Int, length(network.nbrs[v])) + @cast A[_, i] := loc_exp[i] - -function projectors(network::PEPSNetwork, vertex::NTuple{2, Int}) - i, j = vertex - neighbours = ((i, j-1), (i-1, j), (i, j+1), (i+1, j)) - projector.(Ref(network), Ref(vertex), neighbours) + 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 .- minimum(en))) +end -# @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)) - -# dim = zeros(Int, length(network.nbrs[v])) -# @cast A[_, i] := loc_exp[i] - -# # for proj ∈ _get_projectors(network, v) -# # @cast A[(c, γ), σ] |= A[c, σ] * proj[σ, γ] -# # dim[j] = size(pv, 2) -# 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 .- minimum(en))) -# end - -function peps_tensor(::Type{T}, peps::PEPSNetwork, i::Int, j::Int, β::Real) where {T <: Number} +function peps_tensor(::Type{T}, peps::PEPSNetwork, i::Int, j::Int) where {T <: Number} # generate tensors from projectors - A = build_tensor(peps, (i, j), β) + A = generate_tensor(peps, peps.map[i, j]) # include energy - h = build_tensor(peps, (i, j-1), (i, j), β) - v = build_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, σ] B end +peps_tensor(peps::PEPSNetwork, i::Int, j::Int) = peps_tensor(Float64, peps, i, j) -peps_tensor(peps::PEPSNetwork, i::Int, j::Int, β::Real) = peps_tensor(Float64, peps, i, j, β) - -function SpinGlassTensors.PEPSRow(::Type{T}, peps::PEPSNetwork, i::Int, β::Real) where {T <: Number} - ψ = PEPSRow(T, peps.ncols) - for j ∈ 1:peps.ncols - ψ[j] = peps_tensor(T, peps, i, j, β) +function SpinGlassTensors.PEPSRow(::Type{T}, peps::PEPSNetwork, i::Int) where {T <: Number} + ψ = PEPSRow(T, peps.j_max) + for j ∈ 1:peps.j_max + ψ[j] = peps_tensor(T, peps, i, j) end ψ end -SpinGlassTensors.PEPSRow(peps::PEPSNetwork, i::Int, β::Real) = PEPSRow(Float64, peps, i, β) +SpinGlassTensors.PEPSRow(peps::PEPSNetwork, i::Int) = PEPSRow(Float64, peps, i) function SpinGlassTensors.MPO(::Type{T}, peps::PEPSNetwork, i::Int, - β::Real, - config::Dict{NTuple{2, Int}, Int} = Dict{NTuple{2, Int}, Int}() -) where {T <: Number} + config::Dict{Int, Int} = Dict{Int, Int}() + ) where {T <: Number} - W = MPO(T, peps.ncols) - R = PEPSRow(T, peps, i, β) + W = MPO(T, peps.j_max) + R = PEPSRow(T, peps, i) for (j, A) ∈ enumerate(R) - v = get(config, peps.vertex_map((i, j)), nothing) + v = get(config, j + peps.j_max * (i - 1), nothing) if v !== nothing @cast B[l, u, r, d] |= A[l, u, r, d, $(v)] else @@ -152,79 +122,81 @@ end SpinGlassTensors.MPO(peps::PEPSNetwork, i::Int, - β::Real, - config::Dict{NTuple{2, Int}, Int} = Dict{NTuple{2, Int}, Int}() -) = MPO(Float64, peps, i, β, config) - -function compress( - ψ::AbstractMPS, - peps::PEPSNetwork; - bond_dim=typemax(Int), - var_tol=1E-8, - sweeps=4 -) - if bond_dimension(ψ) < bond_dim return ψ end - compress(ψ, bond_dim, var_tol, sweeps) -end + config::Dict{Int, Int} = Dict{Int, Int}() + ) = MPO(Float64, peps, i, config) +function compress(ψ::AbstractMPS, peps::PEPSNetwork) + Dcut = peps.args["bond_dim"] + if bond_dimension(ψ) < Dcut return ψ end + compress(ψ, Dcut, peps.args["var_tol"], peps.args["sweeps"]) +end @memoize function SpinGlassTensors.MPS( peps::PEPSNetwork, i::Int, - β::Real, - cfg::Dict{NTuple{2, Int}, Int} = Dict{NTuple{2, Int}, Int}(); - bond_dim=typemax(Int), - var_tol=1E-8, - sweeps=4 -) - if i > peps.nrows return IdentityMPS() end - W = MPO(peps, i, β, cfg) - ψ = MPS(peps, i+1, β, cfg) - compress(W * ψ, peps, bond_dim=bond_dim, var_tol=var_tol, sweeps=sweeps) + cfg::Dict{Int, Int} = Dict{Int, Int}(), + ) + if i > peps.i_max return IdentityMPS() end + W = MPO(peps, i, cfg) + ψ = MPS(peps, i+1, cfg) + compress(W * ψ, peps) end function contract_network( peps::PEPSNetwork, - β::Real, - config::Dict{NTuple{2, Int}, Int} = Dict{NTuple{2, Int}, Int}(), + config::Dict{Int, Int} = Dict{Int, Int}(), ) - ψ = MPS(peps, 1, β, config) + ψ = MPS(peps, 1, config) prod(dropindices(ψ))[] end - @inline function get_coordinates(peps::PEPSNetwork, k::Int) - ceil(Int, k / peps.ncols), (k - 1) % peps.ncols + 1 + ceil(Int, k / peps.j_max), (k - 1) % peps.j_max + 1 end -function generate_boundary(network::PEPSNetwork, v::NTuple{2, Int}, w::NTuple{2, Int}, state::Int) - if v ∉ vertices(network.network_graph) return 1 end - loc_dim = length(local_energy(network, v)) - pv = projector(network, v, w) +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}, w::NTuple{2, Int}) i, j = w - ∂v = zeros(Int, peps.ncols + 1) + ∂v = zeros(Int, peps.j_max + 1) # on the left below for k ∈ 1:j-1 - ∂v[k] = generate_boundary(peps, (i, k), (i+1, k), _get_local_state(peps, v, (i, k))) + ∂v[k] = generate_boundary( + 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, (i, j-1), (i, j), _get_local_state(peps, v, (i, j-1))) + ∂v[j] = generate_boundary( + 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.ncols - ∂v[k+1] = generate_boundary(peps, (i-1, k), (i, k), _get_local_state(peps, v, (i-1, k))) + 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)) + ) end ∂v end function _get_local_state(peps::PEPSNetwork, σ::Vector{Int}, w::NTuple{2, Int}) - k = w[2] + peps.ncols * (w[1] - 1) + k = w[2] + peps.j_max * (w[1] - 1) 0 < k <= length(σ) ? σ[k] : 1 end @@ -234,16 +206,19 @@ function _normalize_probability(prob::Vector{T}) where {T <: Number} prob / sum(prob) end -function conditional_probability(peps::PEPSNetwork, v::Vector{Int}, β::Real) +function conditional_probability( + peps::PEPSNetwork, + v::Vector{Int}, + ) i, j = get_coordinates(peps, length(v)+1) ∂v = generate_boundary(peps, v, (i, j)) - W = MPO(peps, i, β) - ψ = MPS(peps, i+1, β) + W = MPO(peps, i) + ψ = MPS(peps, i+1) L = left_env(ψ, ∂v[1:j-1]) - R = right_env(ψ, W, ∂v[j+2:peps.ncols+1]) - A = peps_tensor(peps, i, j, β) + R = right_env(ψ, W, ∂v[j+2:peps.j_max+1]) + A = peps_tensor(peps, i, j) l, u = ∂v[j:j+1] M = ψ[j] @@ -253,29 +228,31 @@ function conditional_probability(peps::PEPSNetwork, v::Vector{Int}, β::Real) _normalize_probability(prob) end -function bond_energy(network, u::NTuple{2, Int}, v::NTuple{2, Int}, σ::Int) - fg_u, fg_v = network.vertex_map(u), network.vertex_map(v) - if has_edge(network.factor_graph, fg_u, fg_v) - pu, en, pv = get_prop.(Ref(network.factor_graph), Ref(fg_u), Ref(fg_v), (:pl, :en, :pr)) +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(network.factor_graph, fg_v, fg_u) - pv, en, pu = get_prop.(Ref(network.factor_graph), Ref(fg_v), Ref(fg_u), (:pl, :en, :pr)) + 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(length(local_energy(network, u))) + energies = zeros(get_prop(fg, u, :loc_dim)) end vec(energies) end -function update_energy(network::AbstractGibbsNetwork, σ::Vector{Int}) +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, (i, j), (i, j-1), σil) + - bond_energy(network, (i, j), (i-1, j), σkj) + - local_energy(network, (i, j)) + 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 diff --git a/src/SpinGlassEngine.jl b/src/SpinGlassEngine.jl index d56a7897..76421940 100644 --- a/src/SpinGlassEngine.jl +++ b/src/SpinGlassEngine.jl @@ -19,8 +19,6 @@ function LinearAlgebra.dot(ψ::AbstractMPS, state::Union{Vector, NTuple}) tr(C) end -include("network_operations.jl") -include("network_interface.jl") include("MPS_search.jl") include("search.jl") include("PEPS.jl") diff --git a/src/attic.jl b/src/attic.jl deleted file mode 100644 index 88b46126..00000000 --- a/src/attic.jl +++ /dev/null @@ -1,54 +0,0 @@ -# ρ needs to be ∈ the right canonical form -function solve(ψ::AbstractMPS, keep::Int) - @assert keep > 0 "Number of states has to be > 0" - T = eltype(ψ) - - keep_extra = keep - pCut = prob = 0. - k = 1 - - if keep < prod(rank(ψ)) - keep_extra += 1 - end - - states = fill([], 1, k) - left_env = _make_left_env(ψ, k) - - for (i, M) ∈ enumerate(ψ) - _, d, b = size(M) - - pdo = zeros(T, k, d) - LL = _make_LL(ψ, b, k, d) - config = zeros(Int, i, k, d) - - for j ∈ 1:k - L = left_env[:, :, j] - - for σ ∈ local_basis(d) - m = idx(σ) - LL[:, :, j, m] = M[:, m, :]' * (L * M[:, m, :]) - pdo[j, m] = tr(LL[:, :, j, m]) - config[:, j, m] = vcat(states[:, j]..., σ) - end - end - - perm = collect(1: k * d) - k = min(k * d, keep_extra) - - if k >= keep_extra - partialsortperm!(perm, vec(pdo), 1:k, rev=true) - prob = vec(pdo)[perm] - pCut < last(prob) ? pCut = last(prob) : () - end - - @cast A[α, β, (l, d)] |= LL[α, β, l, d] - left_env = A[:, :, perm] - - @cast B[α, (l, d)] |= config[α, l, d] - states = B[:, perm] - end - states[:, 1:keep], prob[1:keep], pCut -end - -_make_LL(ψ::AbstractMPS, b::Int, k::Int, d::Int) = zeros(eltype(ψ), b, b, k, d) -_make_left_env(ψ::AbstractMPS, k::Int) = ones(eltype(ψ), 1, 1, k) diff --git a/src/network_interface.jl b/src/network_interface.jl deleted file mode 100644 index c8ab5d4c..00000000 --- a/src/network_interface.jl +++ /dev/null @@ -1,97 +0,0 @@ -using LabelledGraphs - -export - AbstractGibbsNetwork, - network_graph, - vertex_map, - projectors, - local_energy, - interaction_energy, - build_tensor - - -# S: type of the vertex of network -# T: type of the vertex of underlying factor graph -abstract type AbstractGibbsNetwork{S, T} end - - -struct NotImplementedError{M} <: Exception - m::M - NotImplementedError(m::M) where {M} = new{M}(m) -end - -not_implmented(name) = throw(NotImplementedError(name)) - - -factor_graph(network::AbstractGibbsNetwork{S, T}) where {S, T} = network.factor_graph - -network_graph(network::AbstractGibbsNetwork{S, T}) where {S, T} = network.network_graph - -vertex_map(network::AbstractGibbsNetwork{S, T}) where {S, T} = network.vertex_map - -projectors(network::AbstractGibbsNetwork{S, T}, vertex::S) where {S, T} = not_implmented("projectors") - - -function projector(network::AbstractGibbsNetwork{S, T}, v::S, w::S) where {S, T} - fg = factor_graph(network) - vmap = vertex_map(network) - fg_v, fg_w = vmap(v), vmap(w) - - if has_edge(fg, fg_w, fg_v) - get_prop(fg, fg_w, fg_v, :pr)' - elseif has_edge(fg, fg_v, fg_w) - get_prop(fg, fg_v, fg_w, :pl) - else - loc_dim = length(local_energy(network, v)) - ones(loc_dim, 1) - end -end - -function spectrum(network::AbstractGibbsNetwork{S, T}, vertex::S) where {S, T} - get_prop(factor_graph(network), vertex_map(network)(vertex), :spectrum) -end - -function local_energy(network::AbstractGibbsNetwork{S, T}, vertex::S) where {S, T} - spectrum(network, vertex).energies -end - -function interaction_energy(network::AbstractGibbsNetwork{S, T}, v::S, w::S) where {S, T} - fg = factor_graph(network) - vmap = vertex_map(network) - fg_v, fg_w = vmap(v), vmap(w) - if has_edge(fg, fg_w, fg_v) - en = get_prop(fg, fg_w, fg_v, :en)' - elseif has_edge(fg, fg_v, fg_w) - en = get_prop(fg, fg_v, fg_w, :en) - else - en = zeros(1, 1) - end - en -end - -@memoize function build_tensor(network::AbstractGibbsNetwork{S, T}, v::S, β::Real) where {S, T} - # TODO: does this require full network, or can we pass only fg? - loc_exp = exp.(-β .* local_energy(network, v)) - - projs = projectors(network, v) - dim = zeros(Int, length(projs)) - @cast A[_, i] := loc_exp[i] - - for (j, pv) ∈ enumerate(projs) - @cast A[(c, γ), σ] |= A[c, σ] * pv[σ, γ] - dim[j] = size(pv, 2) - end - reshape(A, dim..., :) -end - -@memoize function build_tensor(network::AbstractGibbsNetwork{S, T}, v::S, w::S, β::Real) where {S, T} - en = interaction_energy(network, v, w) - exp.(-β .* (en .- minimum(en))) -end - -function is_compatible(factor_graph, network_graph, vertex_map) - all( - has_edge(network_graph, src(edge), dst(edge)) - for edge ∈ edges(factor_graph) - ) -end diff --git a/src/network_operations.jl b/src/network_operations.jl deleted file mode 100644 index 81004282..00000000 --- a/src/network_operations.jl +++ /dev/null @@ -1,84 +0,0 @@ -export LatticeTransformation, rotation, reflection - - -struct LatticeTransformation - permutation::NTuple{4, Int} - flips_dimensions::Bool -end - -Base.:(∘)(op1::LatticeTransformation, op2::LatticeTransformation) = - LatticeTransformation( - op1.permutation[collect(op2.permutation)], - op1.flips_dimensions ⊻ op2.flips_dimensions - ) - - -function reflection(axis::Symbol) - if axis == :x - perm = (4, 3, 2, 1) - flips = false - elseif axis == :y - perm = (2, 1, 4, 3) - flips = false - elseif axis == :diag - perm = (1, 4, 3, 2) - flips = true - elseif axis == :antydiag - perm = (3, 2, 1, 4) - flips = true - else - throw(ArgumentError("Unknown reflection axis: $(axis)")) - end - LatticeTransformation(perm, flips) -end - - -function rotation(θ::Int) - if θ % 90 != 0 - ArgumentError("Only integral multiplicities of 90° can be passed as θ.") - end - θ = θ % 360 - if θ == 0 - LatticeTransformation((1, 2, 3, 4), false) - elseif θ == 90 - LatticeTransformation((4, 1, 2, 3), true) - else - rotation(θ - 90) ∘ rotation(90) - end -end - - -function check_bounds(m, n) - function _check(i, j) - if i < 1 || i > m || j < 1 || j > n - throw(ArgumentError("Point ($(i), $(j)) not in $(m)x$(n) lattice.")) - end - true - end -end - - -function vertex_map(vert_permutation::NTuple{4, Int}, nrows, ncols) - if vert_permutation == (1, 2, 3, 4) # - f = (i, j) -> (i, j) - elseif vert_permutation == (4, 1, 2, 3) # 90 deg rotation - f = (i, j) -> (nrows - j + 1, i) - elseif vert_permutation == (3, 4, 1, 2) # 180 deg rotation - f = (i, j) -> (nrows - i + 1, ncols - j + 1) - elseif vert_permutation == (2, 3, 4, 1) # 270 deg rotation - f = (i, j) -> (j, ncols - i + 1) - elseif vert_permutation == (2, 1, 4, 3) # :y reflection - f = (i, j) -> (i, ncols - j + 1) - elseif vert_permutation == (4, 3, 2, 1) # :x reflection - f = (i, j) -> (nrows - i + 1, j) - elseif vert_permutation == (1, 4, 3, 2) # :diag reflection - f = (i, j) -> (j, i) - elseif vert_permutation == (3, 2, 1, 4) # :antydiag reflection - f = (i, j) -> (nrows -j + 1, ncols - i + 1) - else - ArgumentError("$(vert_permutation) does not define square isometry.") - end - (tuple) -> f(tuple[1], tuple[2]) -end - -vertex_map(trans::LatticeTransformation, m::Int, n::Int) = vertex_map(trans.permutation, m, n) diff --git a/src/search.jl b/src/search.jl index 62eb6c00..13ce3063 100644 --- a/src/search.jl +++ b/src/search.jl @@ -2,7 +2,7 @@ export AbstractGibbsNetwork export low_energy_spectrum export Solution -# abstract type AbstractGibbsNetwork end +abstract type AbstractGibbsNetwork end struct Solution energies::Vector{Float64} @@ -27,9 +27,9 @@ function _bound(probabilities::Vector{Float64}, cut::Int) k = length(probabilities) second_phase = false - if k > cut + 1 + if k > cut + 1 k = cut + 1 - second_phase = true + second_phase = true end idx = partialsortperm(probabilities, 1:k, rev=true) @@ -44,18 +44,17 @@ end function _branch_and_bound( sol::Solution, network::AbstractGibbsNetwork, - node::NTuple{2, Int}, + node::Int, cut::Int, - β::Real ) # branch pdo, eng, cfg = Float64[], Float64[], Vector{Int}[] - k = length(local_energy(network, node)) + k = get_prop(network.fg, node, :loc_dim) for (p, σ, e) ∈ zip(sol.probabilities, sol.states, sol.energies) - pdo = [pdo; p .* conditional_probability(network, σ, β)] + pdo = [pdo; p .* conditional_probability(network, σ)] eng = [eng; e .+ update_energy(network, σ)] cfg = _branch_state(cfg, σ, collect(1:k)) end @@ -71,18 +70,17 @@ end #TODO: incorporate "going back" move to improve alghoritm function low_energy_spectrum( network::AbstractGibbsNetwork, - cut::Int, - β::Real + cut::Int ) sol = Solution([0.], [[]], [1.], -Inf) - perm = zeros(Int, nv(network.factor_graph)) # TODO: to be removed + perm = zeros(Int, nv(network.fg)) # TODO: to be removed #TODO: this should be replaced with the iteration over fg that is consistent with the order network - for i ∈ 1:network.nrows, j ∈ 1:network.ncols - v_fg = network.factor_graph.reverse_label_map[network.vertex_map((i, j))] - perm[v_fg] = j + network.ncols * (i - 1) - sol = _branch_and_bound(sol, network, (i, j), cut, β) + for i ∈ 1:network.i_max, j ∈ 1:network.j_max + v_fg = network.map[i, j] + perm[v_fg] = j + network.j_max * (i - 1) + sol = _branch_and_bound(sol, network, v_fg, cut) end K = partialsortperm(sol.energies, 1:length(sol.energies), rev=false) diff --git a/test/PEPS.jl b/test/PEPS.jl index c28b986a..b022f4c9 100644 --- a/test/PEPS.jl +++ b/test/PEPS.jl @@ -25,7 +25,7 @@ m = 3 n = 4 t = 3 -β = 1.0 +β = 1 L = m * n * t T = Float64 @@ -43,14 +43,13 @@ fg = factor_graph( x, y = m, n -for transform ∈ [rotation.([0, 90, 180, 270])..., reflection.([:x, :y, :diag, :antydiag])...] - peps = PEPSNetwork(x, y, fg, transform) +for origin ∈ (:NW, :SW, :WS, :WN, :NE, :EN, :SE, :ES) + peps = PEPSNetwork(x, y, fg, β, origin) ψ = IdentityMPS() - - for i ∈ peps.nrows:-1:1 - ψ = MPO(T, peps, i, β) * ψ - @test MPS(peps, i, β) ≈ ψ + for i ∈ peps.i_max:-1:1 + ψ = MPO(T, peps, i) * ψ + @test MPS(peps, i) ≈ ψ end end diff --git a/test/contract.jl b/test/contract.jl index 421aad78..0f4089c0 100644 --- a/test/contract.jl +++ b/test/contract.jl @@ -1,8 +1,3 @@ -using LightGraphs -using LabelledGraphs -using MetaGraphs - - @testset "peps_contract correctly collapse the peps network" begin # Grid @@ -26,21 +21,21 @@ using MetaGraphs fg = factor_graph( ig, - Dict((1, 1) => 4, (1, 2) => 2), + Dict(1 => 4, 2 => 2), spectrum = full_spectrum, - cluster_assignment_rule = Dict(1 => (1, 1), 2 => (1, 1), 3 => (1, 2), 4 => (1, 2)) + cluster_assignment_rule = Dict(1 => 1, 2 => 1, 3 => 2, 4 => 2) ) - e, p = get_prop(fg, (1, 1), (1, 2), :en), get_prop(fg, (1, 1), (1, 2), :pr) + e, p = get_prop(fg, 1, 2, :en), get_prop(fg, 1, 2, :pr) ϕ = exp(β * minimum(e * p)) for i ∈ 1:4, j ∈ 1:2 - cfg = Dict((1, 1) => i, (1, 2) => j) + cfg = Dict(1 => i, 2 => j) Z = [] - for transform ∈ [rotation.([0, 90, 180, 270])..., reflection.([:x, :y, :diag, :antydiag])...] - peps = PEPSNetwork(m, n, fg, transform) - p = contract_network(peps, β, cfg) + for origin ∈ (:NW, :SW, :WS, :WN) + peps = PEPSNetwork(m, n, fg, β, origin) + p = contract_network(peps, cfg) push!(Z, p) end @@ -53,6 +48,6 @@ using MetaGraphs ϱ = reshape(ρ, (4, 2)) * ϕ # probabilities should agree - @test first(Z) ≈ ϱ[cfg[(1, 1)], cfg[(1, 2)]] + @test first(Z) ≈ ϱ[cfg[1], cfg[2]] end end diff --git a/test/network_operations.jl b/test/network_operations.jl deleted file mode 100644 index 4d8855b9..00000000 --- a/test/network_operations.jl +++ /dev/null @@ -1,126 +0,0 @@ -using Test - - -@testset "0 degree rotation is and identity" begin - op = vertex_map(rotation(0), 10, 5) - - @test op((1, 1)) == (1, 1) - @test op((10, 5)) == (10, 5) - @test op((1, 5)) == (1, 5) - @test op((10, 1)) == (10, 1) - @test op((4, 3)) == (4, 3) - @test op((2, 2)) == (2, 2) -end - - -@testset "90 degree rotation of square lattice correctly transforms vertices" begin - op = vertex_map(rotation(90), 10, 5) - - @test op((1, 1)) == (10, 1) - @test op((5, 10)) == (1, 5) - @test op((1, 10)) == (1, 1) - @test op((5, 1)) == (10, 5) - @test op((4, 3)) == (8, 4) - @test op((3, 4)) == (7, 3) - @test op((2, 2)) == (9, 2) -end - - -@testset "180 degree rotation of square lattice correctly transforms vertices" begin - op = vertex_map(rotation(180), 10, 5) - - @test op((1, 1)) == (10, 5) - @test op((10, 5)) == (1, 1) - @test op((1, 5)) == (10, 1) - @test op((10, 1)) == (1, 5) - @test op((4, 3)) == (7, 3) -end - - -@testset "270 degree rotation of square lattice correctly transforms vertices" begin - op = vertex_map(rotation(270), 10, 5) - - @test op((1, 1)) == (1, 5) - @test op((5, 1)) == (1, 1) - @test op((5, 10)) == (10, 1) - @test op((1, 10)) == (10, 5) - @test op((4, 3)) == (3, 2) - @test op((3, 4)) == (4, 3) -end - - -@testset "Reflection around x axis correctly transforms vertices" begin - op = vertex_map(reflection(:x), 10, 5) - - @test op((1, 1)) == (10, 1) - @test op((10, 5)) == (1, 5) - @test op((1, 5)) == (10, 5) - @test op((10, 1)) == (1, 1) - @test op((4, 3)) == (7, 3) - @test op((2, 2)) == (9, 2) -end - - -@testset "Reflection around y axis correctly transforms vertices" begin - op = vertex_map(reflection(:y), 10, 5) - - @test op((1, 1)) == (1, 5) - @test op((10, 5)) == (10, 1) - @test op((1, 5)) == (1, 1) - @test op((10, 1)) == (10, 5) - @test op((4, 3)) == (4, 3) - @test op((2, 2)) == (2, 4) -end - - -@testset "Reflection around diag correctly transforms vertices" begin - op = vertex_map(reflection(:diag), 10, 5) - - @test op((1, 1)) == (1, 1) - @test op((5, 10)) == (10, 5) - @test op((1, 10)) == (10, 1) - @test op((5, 1)) == (1, 5) - @test op((4, 3)) == (3, 4) - @test op((2, 2)) == (2, 2) -end - - -@testset "Reflection around antydiag correctly transforms vertices" begin - op = vertex_map(reflection(:antydiag), 10, 5) - - @test op((1, 1)) == (10, 5) - @test op((5, 10)) == (1, 1) - @test op((1, 10)) == (1, 5) - @test op((5, 1)) == (10, 1) - @test op((4, 3)) == (8, 2) - @test op((3, 4)) == (7, 3) - @test op((2, 2)) == (9, 4) -end - - -for transform ∈ (rotation.([0, 90, 180, 270])..., reflection.([:x, :y, :diag, :antydiag])...) -@testset "$(transform) of square lattice is a bijection" begin - op = vertex_map(transform, 3, 3) - all_points = [(i, j) for i ∈ 1:3, j ∈ 1:3] - @test Set(op.(all_points)) == Set(all_points) -end -end - - -for transform ∈ (rotation.([90, 270])..., reflection.([:diag, :antydiag])...) -@testset "$(transform) of rectangular lattice is bijection flipping lattice dimensions" begin - op = vertex_map(transform, 2, 3) - original_points = [(i, j) for i ∈ 1:3, j ∈ 1:2] - transformed_points = [(i, j) for i ∈ 1:2, j ∈ 1:3] - @test Set(op.(original_points)) == Set(transformed_points) -end -end - - -for transform ∈ (rotation.([0, 180])..., reflection.([:x, :y])...) -@testset "$(transform) of rectangular lattice is bijection leaving lattice dimensions unchanged" begin - op = vertex_map(transform, 7, 3) - all_points = [(i, j) for i ∈ 1:7, j ∈ 1:3] - @test Set(op.(all_points)) == Set(all_points) -end -end diff --git a/test/runtests.jl b/test/runtests.jl index 4e4efeb2..2f423706 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -45,8 +45,7 @@ push!(my_tests, "PEPS.jl", "contract.jl", "search.jl", - "search_2.jl", - "network_operations.jl" + "search_2.jl" ) for my_test in my_tests diff --git a/test/search.jl b/test/search.jl index 72ad9461..bc196d29 100644 --- a/test/search.jl +++ b/test/search.jl @@ -1,5 +1,3 @@ -using MetaGraphs - @testset "Simplest possible system of two spins" begin # # ----------------- Ising model ------------------ @@ -39,13 +37,11 @@ using MetaGraphs # construct factor graph with no approx fg = factor_graph( ig, - Dict((1, 1) => 2, (1, 2) => 2), + Dict(1 => 2, 2 => 2), spectrum = full_spectrum, - cluster_assignment_rule = Dict(1 => (1, 1), 2 => (1, 2)), # treat it as a grid with 1 spin cells + cluster_assignment_rule = Dict(1 => 1, 2 => 2), # treat it as a grid with 1 spin cells ) - - # set parameters to contract exactely control_params = Dict( "bond_dim" => typemax(Int), @@ -58,7 +54,7 @@ using MetaGraphs ϱ = gibbs_tensor(ig, β) # split on the bond - p1, e, p2 = get_prop.(Ref(fg), Ref((1, 1)), Ref((1, 2)), (:pl, :en, :pr)) + p1, e, p2 = get_prop.(Ref(fg), 1, 2, (:pl, :en, :pr)) @testset "has correct energy on the bond" begin en = [ J12 * σ * η for σ ∈ [-1, 1], η ∈ [-1, 1]] @@ -66,18 +62,16 @@ using MetaGraphs @test p1 ≈ p2 ≈ I end - for transform ∈ [rotation.([0, 90, 180, 270])..., reflection.([:x, :y, :diag, :antydiag])...] - peps = PEPSNetwork(m, n, fg, transform) - cluster_to_spin = Dict((1, 1) => 1, (1, 2) => 2) - #cluster_to_spin = Dict((1, 1) => 1,(1, 2) => 2) + for origin ∈ (:NW, :SW, :WS, :WN, :NE, :EN, :SE, :ES) + peps = PEPSNetwork(m, n, fg, β, origin, control_params) - @testset "has properly built PEPS tensors given transformation $(transform)" begin + @testset "has properly built PEPS tensors given origin at $(origin)" begin # horizontal alignment - 1 row, 2 columns - if peps.nrows == 1 && peps.ncols == 2 - @test !transform.flips_dimensions + if peps.i_max == 1 && peps.j_max == 2 + @test origin ∈ (:NW, :SW, :SE, :NE) - l, k = cluster_to_spin[peps.vertex_map((1, 1))], cluster_to_spin[peps.vertex_map((1, 2))] + l, k = peps.map[1, 1], peps.map[1, 2] v1 = [exp(-β * D[l, l] * σ) for σ ∈ [-1, 1]] v2 = [exp(-β * D[k, k] * σ) for σ ∈ [-1, 1]] @@ -89,14 +83,14 @@ using MetaGraphs @reduce ρ[σ, η] := sum(l) A[1, 1, l, 1, σ] * B[l, 1, 1, 1, η] if l == 2 ρ = ρ' end - R = PEPSRow(peps, 1, β) + R = PEPSRow(peps, 1) @test [R[1], R[2]] ≈ [A, B] # vertical alignment - 1 column, 2 rows - elseif peps.nrows == 2 && peps.ncols == 1 - @test transform.flips_dimensions - l, k = cluster_to_spin[peps.vertex_map((1, 1))], cluster_to_spin[peps.vertex_map((2, 1))] - # l, k = peps.map[1, 1], peps.map[2, 1] + elseif peps.i_max == 2 && peps.j_max == 1 + @test origin ∈ (:WN, :WS, :ES, :EN) + + l, k = peps.map[1, 1], peps.map[2, 1] v1 = [exp(-β * D[l, l] * σ) for σ ∈ [-1, 1]] v2 = [exp(-β * D[k, k] * σ) for σ ∈ [-1, 1]] @@ -108,8 +102,8 @@ using MetaGraphs @reduce ρ[σ, η] := sum(u) A[1, 1, 1, u, σ] * B[1, u, 1, 1, η] if l == 2 ρ = ρ' end - @test PEPSRow(peps, 1, β)[1] ≈ A - @test PEPSRow(peps, 2, β)[1] ≈ B + @test PEPSRow(peps, 1)[1] ≈ A + @test PEPSRow(peps, 2)[1] ≈ B end @testset "which produces correct Gibbs state" begin @@ -118,12 +112,12 @@ using MetaGraphs end # solve the problem using B & B - sol = low_energy_spectrum(peps, num_states, β) + sol = low_energy_spectrum(peps, num_states) - @testset "has correct spectrum given transformation $(transform)" begin + @testset "has correct spectrum given the origin at $(origin)" begin for (σ, η) ∈ zip(exact_spectrum.states, sol.states) - for i ∈ 1:peps.nrows, j ∈ 1:peps.ncols - v = j + peps.ncols * (i - 1) + for i ∈ 1:peps.i_max, j ∈ 1:peps.j_max + v = j + peps.j_max * (i - 1) # 1 --> -1 and 2 --> 1 @test (η[v] == 1 ? -1 : 1) == σ[v] end diff --git a/test/search_2.jl b/test/search_2.jl index 7f74d9de..22861408 100644 --- a/test/search_2.jl +++ b/test/search_2.jl @@ -58,11 +58,11 @@ 21 => 5, 22 => 5, ) - # control_params = Dict( - # "bond_dim" => typemax(Int), - # "var_tol" => 1E-8, - # "sweeps" => 4. - # ) + control_params = Dict( + "bond_dim" => typemax(Int), + "var_tol" => 1E-8, + "sweeps" => 4. + ) instance = "$(@__DIR__)/instances/pathological/test_$(m)_$(n)_$(t).txt" @@ -74,13 +74,13 @@ cluster_assignment_rule=super_square_lattice((m, n, t)) ) - for transform ∈ [rotation.([0, 90, 180, 270])..., reflection.([:x, :y, :diag, :antydiag])...] - peps = PEPSNetwork(m, n, fg, transform) + for origin ∈ (:NW, :SW, :WS, :WN, :NE, :EN, :SE, :ES) + peps = PEPSNetwork(m, n, fg, β, origin, control_params) # solve the problem using B & B - sol = low_energy_spectrum(peps, num_states, β) + sol = low_energy_spectrum(peps, num_states) - @testset "has correct spectrum given the transformation $(transform)" begin + @testset "has correct spectrum given the origin at $(origin)" begin @test sol.energies ≈ exact_energies for (i, σ) ∈ enumerate(sol.states) @test σ ∈ exact_states[deg[i]]