Remove superflous data structures in network and PEPS.jl

src/MPS_search.jl

@@ -10,26 +10,6 @@ struct MPSControl
     dβ::Number
 end
 
-#= it will be used instead of MPSControl
-function _set_control_parameters(
-    args_override::Dict{String, Number}=Dict{String, Number}()
-    )
-    args = Dict(
-        "max_bond" => typemax(Int),
-        "var_ϵ" => 1E-8,
-        "max_sweeps" => 4.,
-        "β" => 1.,
-        "dβ" => 0.01
-    )
-    for k in keys(args_override)
-        str = get(args_override, k, nothing)
-        if str !== nothing push!(args, str) end
-    end
-    args
-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)

src/PEPS.jl

@@ -1,29 +1,25 @@
-export PepsNetwork, contract_network
+export PEPSNetwork, contract_network
 export MPO, MPS, generate_boundary
 
-function _set_control_parameters(
-    args_override::Dict{String, Number}=Dict{String, Number}()
-    )
-    # put here more parameters if needs be
-    args = Dict(
-        "bond_dim" => typemax(Int),
-        "var_tol" => 1E-8,
-        "sweeps" => 4.,
-        "β" => 1.
-    )
-    merge(args, args_override)
-end
+const DEFAULT_CONTROL_PARAMS = Dict(
+    "bond_dim" => typemax(Int),
+    "var_tol" => 1E-8,
+    "sweeps" => 4.,
+    "β" => 1.
+)
 
-mutable struct PepsNetwork
+mutable struct PEPSNetwork
     size::NTuple{2, Int}
     map::Dict
-    network_graph::NetworkGraph
+    fg::MetaDiGraph
+    nbrs::Dict
     origin::Symbol
     i_max::Int
     j_max::Int
+    β::Number
     args::Dict{String, Number}
 
-    function PepsNetwork(
+    function PEPSNetwork(
         m::Int,
         n::Int,
         fg::MetaDiGraph,
@@ -42,28 +38,30 @@ mutable struct PepsNetwork
             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.network_graph = NetworkGraph(fg, nbrs, β)
-        pn.args = _set_control_parameters(args_override)
+        pn.fg = fg
+        pn.nbrs = nbrs
+        pn.β = β
+        pn.args = merge(DEFAULT_CONTROL_PARAMS, args_override)
         pn
     end
 end
 
-generate_tensor(pn::PepsNetwork,
+generate_tensor(pn::PEPSNetwork,
                 m::NTuple{2,Int},
-                ) = generate_tensor(pn.network_graph, pn.map[m])
+                ) = generate_tensor(pn, pn.map[m])
 
-generate_tensor(pn::PepsNetwork,
+generate_tensor(pn::PEPSNetwork,
                 m::NTuple{2, Int},
                 n::NTuple{2, Int},
-                ) = generate_tensor(pn.network_graph, pn.map[m], pn.map[n])
+                ) = generate_tensor(pn, pn.map[m], pn.map[n])
 
-generate_boundary(pn::PepsNetwork,
+generate_boundary(pn::PEPSNetwork,
                   m::NTuple{2, Int},
                   n::NTuple{2, Int},
                   σ::Int,
-                 ) = generate_boundary(pn.network_graph, pn.map[m], pn.map[n], σ)
+                 ) = generate_boundary(pn, pn.map[m], pn.map[n], σ)
 
-function PEPSRow(::Type{T}, peps::PepsNetwork, i::Int) where {T <: Number}
+function PEPSRow(::Type{T}, peps::PEPSNetwork, i::Int) where {T <: Number}
     ψ = PEPSRow(T, peps.j_max)
 
     # generate tensors from projectors
@@ -81,10 +79,10 @@ function PEPSRow(::Type{T}, peps::PepsNetwork, i::Int) where {T <: Number}
     end
     ψ
 end
-PEPSRow(peps::PepsNetwork, i::Int) = PEPSRow(Float64, peps, i)
+PEPSRow(peps::PEPSNetwork, i::Int) = PEPSRow(Float64, peps, i)
 
 function MPO(::Type{T},
-    peps::PepsNetwork,
+    peps::PEPSNetwork,
     i::Int,
     config::Dict{Int, Int} = Dict{Int, Int}()
     ) where {T <: Number}
@@ -103,19 +101,19 @@ function MPO(::Type{T},
     end
     W
 end
-MPO(peps::PepsNetwork,
+MPO(peps::PEPSNetwork,
     i::Int,
     config::Dict{Int, Int} = Dict{Int, Int}()
     ) = MPO(Float64, peps, i, config)
 
-function compress(ψ::AbstractMPS, peps::PepsNetwork)
+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 MPS(
-    peps::PepsNetwork,
+    peps::PEPSNetwork,
     i::Int,
     cfg::Dict{Int, Int} = Dict{Int, Int}(),
     )
@@ -126,22 +124,22 @@ end
 end
 
 function contract_network(
-    peps::PepsNetwork,
+    peps::PEPSNetwork,
     config::Dict{Int, Int} = Dict{Int, Int}(),
     )
     ψ = MPS(peps, 1, config)
     prod(dropdims(ψ))[]
 end
 
 @inline function _get_coordinates(
-    peps::PepsNetwork,
+    peps::PEPSNetwork,
     k::Int
     )
     ceil(k / peps.j_max), (k - 1) % peps.j_max + 1
 end
 
 @inline function _get_local_state(
-    peps::PepsNetwork,
+    peps::PEPSNetwork,
     v::Vector{Int},
     i::Int,
     j::Int,
@@ -152,7 +150,7 @@ end
 end
 
 function generate_boundary(
-    peps::PepsNetwork,
+    peps::PEPSNetwork,
     v::Vector{Int},
     i::Int,
     j::Int,
@@ -187,7 +185,7 @@ function generate_boundary(
 end
 
 function generate_boundary(
-    peps::PepsNetwork,
+    peps::PEPSNetwork,
     v::Vector{Int},
     )
     i, j = _get_coordinates(peps, length(v)+1)
@@ -216,7 +214,7 @@ function _normalize_probability(prob::Vector{T}) where {T <: Number}
 end
 
 function conditional_probability(
-    peps::PepsNetwork,
+    peps::PEPSNetwork,
     v::Vector{Int},
     )
     i, j = _get_coordinates(peps, length(v)+1)

src/SpinGlassPEPS.jl

@@ -22,8 +22,8 @@ module SpinGlassPEPS
     include("ising.jl")
     include("exact.jl")
     include("factor.jl")
-    include("network.jl")
     include("PEPS.jl")
+    include("network.jl")
     include("MPS_search.jl")
     include("notation.jl")
     include("peps_no_types.jl")

src/factor.jl

@@ -4,14 +4,28 @@ export factor_graph, rank_reveal, projectors, split_into_clusters
 function split_into_clusters(vertices, assignment_rule)
     # TODO: check how to do this in functional-style
     clusters = Dict(
-        i => Set{Int}() for i in values(assignment_rule)
+        i => [] for i in values(assignment_rule)
     )
     for v in vertices
         push!(clusters[assignment_rule[v]], v)
     end
     clusters
 end
 
+function split_into_clusters(ig::MetaGraph, assignment_rule)
+    cluster_id_to_verts = Dict(
+        i => Int[] for i in values(assignment_rule)
+    )
+
+    for (i, v) in enumerate(nodes(ig))
+        push!(cluster_id_to_verts[assignment_rule[v]], i)
+    end
+
+    return Dict(
+        i => cluster(ig, verts) for (i, verts) ∈ cluster_id_to_verts
+    )
+end
+
 function factor_graph(
     ig::MetaGraph,
     num_states_cl::Int;
@@ -37,51 +51,39 @@ function factor_graph(
     cluster_assignment_rule::Dict{Int, Int} # e.g. square lattice
 )
     L = maximum(values(cluster_assignment_rule))
+    fg = MetaDiGraph(L)
 
-    fg = MetaDiGraph(L, 0.0)
-
-    cluster_to_verts = split_into_clusters(vertices(ig), cluster_assignment_rule)
-
-    for (v, verts) ∈ cluster_to_verts
-        cl = cluster(ig, verts)
+    for (v, cl) ∈ split_into_clusters(ig, cluster_assignment_rule)
         set_prop!(fg, v, :cluster, cl)
-        r = prod(rank_vec(cl))
-        num_states = get(num_states_cl, v, r)
-        sp = spectrum(cl, num_states=num_states)
+        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))
     end
 
     for i ∈ 1:L, j ∈ i+1:L
-        v = get_prop(fg, i, :cluster)
-        w = get_prop(fg, j, :cluster)
+        v, w = get_prop(fg, i, :cluster), get_prop(fg, j, :cluster)
 
         outer_edges, J = inter_cluster_edges(ig, v, w)
 
         if !isempty(outer_edges)
-            e = SimpleEdge(i, j)
-
-            add_edge!(fg, e)
-            set_prop!(fg, e, :outer_edges, outer_edges)
-
-            st_v = get_prop(fg, i, :spectrum).states
-            st_w = get_prop(fg, j, :spectrum).states
-
-            en = zeros(length(st_v), length(st_w))
-
-            for (k, η) ∈ enumerate(vec(st_w))
-                en[:, k] = energy.(vec(st_v), Ref(J), Ref(η))
-            end
+            en = inter_cluster_energy(
+                get_prop(fg, i, :spectrum).states, J, get_prop(fg, j, :spectrum).states
+            )
 
             pl, en = rank_reveal(en, :PE)
             en, pr = rank_reveal(en, :EP)
 
-            set_prop!(fg, e, :split, (pl, en, pr))
+            add_edge!(
+                fg, i, j,
+                Dict(:outer_edges => outer_edges, :pl => pl, :en => en, :pr => pr)
+            )
         end
     end
     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)

src/ising.jl

@@ -1,6 +1,6 @@
 export ising_graph
 export energy, rank_vec
-export Spectrum, cluster, rank
+export Spectrum, cluster, rank, nodes, basis_size
 
 const Instance = Union{String, Dict}
 const SimpleEdge = LightGraphs.SimpleGraphs.SimpleEdge
@@ -11,7 +11,7 @@ struct Spectrum
     states::Array{Vector{<:Number}}
 end
 
-unique_vertices(ising_tuples) = sort(collect(Set(Iterators.flatten((i, j) for (i, j, _) ∈ ising_tuples))))
+unique_nodes(ising_tuples) = sort(collect(Set(Iterators.flatten((i, j) for (i, j, _) ∈ ising_tuples))))
 
 """
 $(TYPEDSIGNATURES)
@@ -34,28 +34,26 @@ function ising_graph(
         ising = [ (i, j, J) for ((i, j), J) ∈ instance ]
     end
 
-    original_verts = unique_vertices(ising)
-    verts_map = Dict(w => i for (i, w) ∈ enumerate(original_verts))
+    original_nodes = unique_nodes(ising)
+    L = length(original_nodes)
+    nodes_to_vertices = Dict(w => i for (i, w) ∈ enumerate(original_nodes))
 
-    L = length(original_verts)
+    ig = MetaGraph(length(original_nodes))
 
-    ig = MetaGraph(L)
-
-    foreach(args -> set_prop!(ig, args[1], :orig_vert, args[2]), enumerate(original_verts))
+    foreach(args -> set_prop!(ig, args[1], :node, args[2]), enumerate(original_nodes))
 
     J = zeros(L, L)
     h = zeros(L)
 
     # setup the model (J_ij, h_i)
     for (_i, _j, v) ∈ ising
-        i, j = verts_map[_i], verts_map[_j]
+        i, j = nodes_to_vertices[_i], nodes_to_vertices[_j]
         v *= sgn
 
         if i == j
             h[i] = v
         else
-            add_edge!(ig, i, j) &&
-            set_prop!(ig, i, j, :J, v) || throw(ArgumentError("Duplicate Egde ($i, $j)"))
+            add_edge!(ig, i, j, :J, v) || throw(ArgumentError("Duplicate Egde ($i, $j)"))
             J[i, j] = v
         end
     end
@@ -66,19 +64,21 @@ function ising_graph(
         ig,
         :rank,
         Dict{Int, Int}(
-            verts_map[v] => get(rank_override, v, 2) for v in original_verts
+            v => get(rank_override, w, 2) for (w, v) in nodes_to_vertices
         )
     )
 
     set_prop!(ig, :J, J)
     set_prop!(ig, :h, h)
-
+    set_prop!(ig, :nodes_map, nodes_to_vertices)
     ig
 end
 
+nodes(ig::MetaGraph) = collect(get_prop.(Ref(ig), vertices(ig), :node))
 rank_vec(ig::MetaGraph) = collect(values(get_prop(ig, :rank)))
+basis_size(ig::MetaGraph) = prod(prod(rank_vec(ig)))
 
-function cluster(ig::MetaGraph, verts::Set{Int})
+function cluster(ig::MetaGraph, verts)
     sub_ig, vmap = induced_subgraph(ig, collect(verts))
 
     h = get_prop.(Ref(sub_ig), vertices(sub_ig), :h)
@@ -120,3 +120,10 @@ E = -\\sum_<i,j> s_i J_{ij} * s_j - \\sum_j h_i s_j.
 energy(σ::Vector, J::Matrix, η::Vector=σ) = dot(σ, J, η)
 energy(σ::Vector, h::Vector) = dot(h, σ)
 energy(σ::Vector, ig::MetaGraph) = energy(σ, get_prop(ig, :J)) + energy(σ, get_prop(ig, :h))
+
+
+# Please don't make the below another energy method.
+# There is already so much mess going on :)
+function inter_cluster_energy(cl1_states, J::Matrix, cl2_states)
+    hcat(collect.(cl1_states)...)' * J * hcat(collect.(cl2_states)...)
+end