From 3de25f700da2bd7da9ff2c95a04756bf333d7e6c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C5=81ukasz=20Pawela?= Date: Wed, 6 Mar 2024 19:15:34 +0100 Subject: [PATCH] reformat --- benchmarks/spectrum.jl | 4 +- docs/make.jl | 13 +- examples/bp.jl | 82 ++++---- examples/temp.jl | 44 +++-- src/SpinGlassNetworks.jl | 42 ++-- src/bp.jl | 84 +++++--- src/clustered_hamiltonian.jl | 145 ++++++++------ src/ising.jl | 41 ++-- src/lattice.jl | 151 +++++++++------ src/spectrum.jl | 33 ++-- src/truncate.jl | 69 ++++--- src/utils.jl | 52 +++-- test/bp_1site.jl | 251 ++++++++++++------------ test/bp_2site.jl | 349 +++++++++++++++++----------------- test/clustered_hamiltonian.jl | 295 +++++++++++++++------------- test/ising.jl | 136 ++++++++----- test/runtests.jl | 6 +- test/spectrum.jl | 4 +- test/utils.jl | 6 +- 19 files changed, 1026 insertions(+), 781 deletions(-) diff --git a/benchmarks/spectrum.jl b/benchmarks/spectrum.jl index cc241da..e08b3a2 100644 --- a/benchmarks/spectrum.jl +++ b/benchmarks/spectrum.jl @@ -1,10 +1,10 @@ using SpinGlassNetworks -function bench(instance::String, size::NTuple{3, Int}, max_states::Int=100) +function bench(instance::String, size::NTuple{3,Int}, max_states::Int = 100) ig = ising_graph(instance) cl = split_into_clusters(ig, super_square_lattice(size)) - @time sp = brute_force(cl[1, 1], num_states=max_states) + @time sp = brute_force(cl[1, 1], num_states = max_states) nothing end diff --git a/docs/make.jl b/docs/make.jl index 21f5b06..ae35025 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -6,19 +6,18 @@ _pages = [ "Lattice geometries" => "lattice.md", "Clustered hamiltonian" => "clh.md", "Local dimensional reduction" => "bp.md", - "API Reference for auxiliary functions" => "api.md" + "API Reference for auxiliary functions" => "api.md", ] # ============================ -format = Documenter.HTML(edit_link = "master", - prettyurls = get(ENV, "CI", nothing) == "true", -) +format = + Documenter.HTML(edit_link = "master", prettyurls = get(ENV, "CI", nothing) == "true") # format = Documenter.LaTeX(platform="none") makedocs( - sitename="SpinGlassNetworks.jl", + sitename = "SpinGlassNetworks.jl", modules = [SpinGlassNetworks], pages = _pages, - format = format - ) \ No newline at end of file + format = format, +) diff --git a/examples/bp.jl b/examples/bp.jl index 53cdca3..b580145 100644 --- a/examples/bp.jl +++ b/examples/bp.jl @@ -14,56 +14,56 @@ Instance below looks like this: 7 -- 8 -- 9 """ function create_larger_example_clustered_hamiltonian_tree() - instance = Dict( - (1, 1) => 0.5, - (2, 2) => 0.25, - (3, 3) => 0.3, - (4, 4) => 0.1, - (5, 5) => -0.1, - (6, 6) => 0.1, - (7, 7) => 0.0, - (8, 8) => 0.1, - (9, 9) => 0.01, - (1, 2) => -1.0, - (2, 3) => 1.0, - (1, 4) => 1.0, - (4, 5) => 1.0, - (5, 6) => 1.0, - (4, 7) => 1.0, - (7, 8) => 1.0, - (8, 9) => 1.0 - ) + instance = Dict( + (1, 1) => 0.5, + (2, 2) => 0.25, + (3, 3) => 0.3, + (4, 4) => 0.1, + (5, 5) => -0.1, + (6, 6) => 0.1, + (7, 7) => 0.0, + (8, 8) => 0.1, + (9, 9) => 0.01, + (1, 2) => -1.0, + (2, 3) => 1.0, + (1, 4) => 1.0, + (4, 5) => 1.0, + (5, 6) => 1.0, + (4, 7) => 1.0, + (7, 8) => 1.0, + (8, 9) => 1.0, + ) - ig = ising_graph(instance) + ig = ising_graph(instance) - assignment_rule = Dict( - 1 => (1, 1, 1), - 2 => (1, 2, 1), - 3 => (1, 3, 1), - 4 => (2, 1, 1), - 5 => (2, 2, 1), - 6 => (2, 3, 1), - 7 => (3, 1, 1), - 8 => (3, 2, 1), - 9 => (3, 3, 1) - ) + assignment_rule = Dict( + 1 => (1, 1, 1), + 2 => (1, 2, 1), + 3 => (1, 3, 1), + 4 => (2, 1, 1), + 5 => (2, 2, 1), + 6 => (2, 3, 1), + 7 => (3, 1, 1), + 8 => (3, 2, 1), + 9 => (3, 3, 1), + ) - cl_h = clustered_hamiltonian( - ig, - Dict{NTuple{3, Int}, Int}(), - spectrum = full_spectrum, - cluster_assignment_rule = assignment_rule, - ) + cl_h = clustered_hamiltonian( + ig, + Dict{NTuple{3,Int},Int}(), + spectrum = full_spectrum, + cluster_assignment_rule = assignment_rule, + ) - ig, cl_h + ig, cl_h end ig, cl_h = create_larger_example_clustered_hamiltonian_tree() beta = 0.1 iter = 0 -beliefs = belief_propagation(cl_h, beta; iter=iter) +beliefs = belief_propagation(cl_h, beta; iter = iter) for v in vertices(cl_h) - en = get_prop(cl_h, v, :spectrum).energies - println("vertex ", v, " energy = ", en .- minimum(en), " bp = ", beliefs[v]) + en = get_prop(cl_h, v, :spectrum).energies + println("vertex ", v, " energy = ", en .- minimum(en), " bp = ", beliefs[v]) end diff --git a/examples/temp.jl b/examples/temp.jl index 18c512a..78d325b 100644 --- a/examples/temp.jl +++ b/examples/temp.jl @@ -26,7 +26,7 @@ function load_openGM(fname::String, Nx::Integer, Ny::Integer) nn = pop!(F) n = [] - for _ in 1:nn + for _ = 1:nn tt = pop!(F) ny, nx = divrem(tt, Nx) push!(n, ny, nx) @@ -61,17 +61,23 @@ function load_openGM(fname::String, Nx::Integer, Ny::Integer) nn = pop!(J) n = [] - for _ in 1:nn + for _ = 1:nn push!(n, pop!(J)) end upper = lower + prod(n) - functions[ii] = reshape(V[lower + 1:upper], reverse(n)...)' + functions[ii] = reshape(V[lower+1:upper], reverse(n)...)' lower = upper end - result = Dict("fun" => functions, "fac" => factors, "N" => reshape(N, (Ny, Nx)), "Nx" => Nx, "Ny" => Ny) + result = Dict( + "fun" => functions, + "fac" => factors, + "N" => reshape(N, (Ny, Nx)), + "Nx" => Nx, + "Ny" => Ny, + ) result end @@ -85,22 +91,39 @@ function clustered_hamiltonian(fname::String, Nx::Integer = 240, Ny::Integer = 3 cl_h = LabelledGraph{MetaDiGraph}(sort(collect(values(clusters)))) for v ∈ cl_h.labels x, y = v - sp = Spectrum(Vector{Real}(undef, 1), Array{Vector{Int}}(undef, 1, 1), Vector{Int}(undef, 1)) + sp = Spectrum( + Vector{Real}(undef, 1), + Array{Vector{Int}}(undef, 1, 1), + Vector{Int}(undef, 1), + ) set_props!(cl_h, v, Dict(:cluster => v, :spectrum => sp)) end for (index, value) in factors if length(index) == 2 y, x = index Eng = sum(functions[value]) - set_props!(cl_h, (x+1, y+1), Dict(:eng => Eng)) + set_props!(cl_h, (x + 1, y + 1), Dict(:eng => Eng)) elseif length(index) == 4 y1, x1, y2, x2 = index add_edge!(cl_h, (x1 + 1, y1 + 1), (x2 + 1, y2 + 1)) - Eng = sum(functions[value], dims=2) - set_props!(cl_h, (x1 + 1, y1 + 1), (x2 + 1, y2 + 1), Dict(:outer_edges=> ((x1 + 1, y1 + 1), (x2 + 1, y2 + 1)), - :eng => Eng, :pl => I, :pr => I)) + Eng = sum(functions[value], dims = 2) + set_props!( + cl_h, + (x1 + 1, y1 + 1), + (x2 + 1, y2 + 1), + Dict( + :outer_edges => ((x1 + 1, y1 + 1), (x2 + 1, y2 + 1)), + :eng => Eng, + :pl => I, + :pr => I, + ), + ) else - throw(ErrorException("Something is wrong with factor index, it has length $(length(index))")) + throw( + ErrorException( + "Something is wrong with factor index, it has length $(length(index))", + ), + ) end end @@ -113,4 +136,3 @@ filename = "/home/tsmierzchalski/.julia/dev/SpinGlassNetworks/examples/penguin-s cf = clustered_hamiltonian(filename, x, y) - diff --git a/src/SpinGlassNetworks.jl b/src/SpinGlassNetworks.jl index aa1dd8d..75e9054 100644 --- a/src/SpinGlassNetworks.jl +++ b/src/SpinGlassNetworks.jl @@ -1,25 +1,25 @@ module SpinGlassNetworks - using LabelledGraphs - using Graphs - using MetaGraphs # TODO: to be replaced by MetaGraphsNext - using CSV - using JLD2 - using DocStringExtensions - using LinearAlgebra, MKL - using Base.Cartesian - using SparseArrays - using HDF5 - using CUDA, CUDA.CUSPARSE - using SpinGlassTensors - import Base.Prehashed +using LabelledGraphs +using Graphs +using MetaGraphs # TODO: to be replaced by MetaGraphsNext +using CSV +using JLD2 +using DocStringExtensions +using LinearAlgebra, MKL +using Base.Cartesian +using SparseArrays +using HDF5 +using CUDA, CUDA.CUSPARSE +using SpinGlassTensors +import Base.Prehashed - include("ising.jl") - include("spectrum.jl") - include("lattice.jl") - #include("projectors.jl") - include("clustered_hamiltonian.jl") - include("bp.jl") - include("truncate.jl") - include("utils.jl") +include("ising.jl") +include("spectrum.jl") +include("lattice.jl") +#include("projectors.jl") +include("clustered_hamiltonian.jl") +include("bp.jl") +include("truncate.jl") +include("utils.jl") end # module diff --git a/src/bp.jl b/src/bp.jl index 92ab0f4..1ef02ef 100644 --- a/src/bp.jl +++ b/src/bp.jl @@ -1,6 +1,5 @@ -export - belief_propagation, +export belief_propagation, clustered_hamiltonian_2site, projector, get_neighbors, @@ -9,7 +8,7 @@ export merge_vertices_cl_h, local_energy, interaction_energy, - SparseCSC + SparseCSC """ $(TYPEDSIGNATURES) @@ -30,7 +29,12 @@ Belief propagation is an iterative algorithm that computes beliefs by passing me The algorithm continues until convergence or until the specified maximum number of iterations is reached. The beliefs are computed based on the inverse temperature parameter `beta`, which controls the influence of energy values on the beliefs. """ -function belief_propagation(cl_h::LabelledGraph{S, T}, beta::Real; tol=1e-6, iter=1) where {S, T} +function belief_propagation( + cl_h::LabelledGraph{S,T}, + beta::Real; + tol = 1e-6, + iter = 1, +) where {S,T} messages_ve = Dict() messages_ev = Dict() @@ -57,7 +61,9 @@ function belief_propagation(cl_h::LabelledGraph{S, T}, beta::Real; tol=1e-6, ite E_local = get_prop(cl_h, v, :spectrum).energies temp = exp.(-(E_local .- minimum(E_local)) * beta) for (n2, pv2, _) in get_neighbors(cl_h, v) - if n1 == n2 continue end + if n1 == n2 + continue + end temp .*= node_messages[n2, v] # messages_ev[n2, v][pv2] end temp ./= sum(temp) @@ -73,7 +79,10 @@ function belief_propagation(cl_h::LabelledGraph{S, T}, beta::Real; tol=1e-6, ite end # Check convergence - converged = all([all(abs.(old_messages_ev[v] .- messages_ev[v]) .< tol) for v in keys(messages_ev)]) + converged = all([ + all(abs.(old_messages_ev[v] .- messages_ev[v]) .< tol) for + v in keys(messages_ev) + ]) end beliefs = Dict() @@ -83,7 +92,7 @@ function belief_propagation(cl_h::LabelledGraph{S, T}, beta::Real; tol=1e-6, ite for (n, pv, _) ∈ get_neighbors(cl_h, v) beliefs[v] .*= messages_ev[n, v][pv] end - beliefs[v] = -log.(beliefs[v])./beta + beliefs[v] = -log.(beliefs[v]) ./ beta beliefs[v] = beliefs[v] .- minimum(beliefs[v]) end @@ -109,7 +118,7 @@ This function retrieves the neighbors of a given vertex in a clustered Hamiltoni It iterates through the edges of the graph and identifies edges connected to the specified vertex. For each neighboring edge, it extracts and returns the neighboring vertex, the associated projector, and the energy. """ -function get_neighbors(cl_h::LabelledGraph{S, T}, vertex::NTuple) where {S, T} +function get_neighbors(cl_h::LabelledGraph{S,T}, vertex::NTuple) where {S,T} neighbors = [] for edge in edges(cl_h) src_node, dst_node = src(edge), dst(edge) @@ -146,7 +155,7 @@ Each field of the `MergedEnergy` struct stores energy values as an `AbstractMatr where `T` is a subtype of the `Real` abstract type. The specific organization and interpretation of these energy values depend on the context in which this struct is used. """ -struct MergedEnergy{T <: Real} +struct MergedEnergy{T<:Real} e11::AbstractMatrix{T} e12::AbstractMatrix{T} e21::AbstractMatrix{T} @@ -220,25 +229,25 @@ function update_message(E_bond::MergedEnergy, message::Vector, beta::Real) R = reshape(reshape(R, sl1 * sr1, sr2) * e22', sl1, sr1, sl2) # [l1, r1, l2] R .*= reshape(e11, sl1, sr1, 1) # [l1, r1, l2] .* [l1, r1, :] R .*= reshape(e21', 1, sr1, sl2) # [l1, r1, l2] .* [:, r1, l2] - R = reshape(sum(R, dims=2), sl1 * sl2) + R = reshape(sum(R, dims = 2), sl1 * sl2) elseif sl1 <= sl2 && sr2 <= sr1 R = reshape(e11', sr1, sl1, 1) .* reshape(message, sr1, 1, sr2) R = reshape(e21 * reshape(R, sr1, sl1 * sr2), sl2, sl1, sr2) R .*= reshape(e12, 1, sl1, sr2) # [l2, l1, r2] .* [:, l1, r2] R .*= reshape(e22, sl2, 1, sr2) # [l2, l1, r2] .* [l2, :, r2] - R = reshape(reshape(sum(R, dims=3), sl2, sl1)', sl1 * sl2) + R = reshape(reshape(sum(R, dims = 3), sl2, sl1)', sl1 * sl2) elseif sl2 <= sl1 && sr1 <= sr2 R = reshape(e22, sl2, 1, sr2) .* reshape(message, 1, sr1, sr2) R = reshape(reshape(R, sl2 * sr1, sr2) * e12', sl2, sr1, sl1) # [l2, r1, l1] R .*= reshape(e11', 1, sr1, sl1) # [l2, r1, l1] .* [:, r1, l1] R .*= reshape(e21, sl2, sr1, 1) # [l2, r1, l1] .* [l2, r1, :] - R = reshape(reshape(sum(R, dims=2), sl2, sl1)', sl1 * sl2) + R = reshape(reshape(sum(R, dims = 2), sl2, sl1)', sl1 * sl2) else # sl2 <= sl1 && sr2 <= sr1 R = reshape(e21', sr1, sl2, 1) .* reshape(message, sr1, 1, sr2) R = reshape(e11 * reshape(R, sr1, sl2 * sr2), sl1, sl2, sr2) R .*= reshape(e12, sl1, 1, sr2) # [l1, l2, r2] .* [l1, :, r2] R .*= reshape(e22, 1, sl2, sr2) # [l1, l2, r2] .* [:, l2, r2] - R = reshape(sum(R, dims=3), sl1 * sl2) + R = reshape(sum(R, dims = 3), sl1 * sl2) end R end @@ -262,7 +271,7 @@ The resulting `new_cl_h` graph represents the 2-site cluster Hamiltonian with si The energy values, projectors, and spectra associated with the new vertices and edges are computed based on the provided temperature parameter `beta`. """ -function clustered_hamiltonian_2site(cl_h::LabelledGraph{S, T}, beta::Real) where {S, T} +function clustered_hamiltonian_2site(cl_h::LabelledGraph{S,T}, beta::Real) where {S,T} unified_vertices = unique([vertex[1:2] for vertex in vertices(cl_h)]) new_cl_h = LabelledGraph{MetaDiGraph}(unified_vertices) @@ -271,7 +280,9 @@ function clustered_hamiltonian_2site(cl_h::LabelledGraph{S, T}, beta::Real) wher vertx = Set() for v in vertices(cl_h) i, j, _ = v - if (i, j) ∈ vertx continue end + if (i, j) ∈ vertx + continue + end E1 = local_energy(cl_h, (i, j, 1)) E2 = local_energy(cl_h, (i, j, 2)) E = energy_2site(cl_h, i, j) .+ reshape(E1, :, 1) .+ reshape(E2, 1, :) @@ -282,12 +293,16 @@ function clustered_hamiltonian_2site(cl_h::LabelledGraph{S, T}, beta::Real) wher edge_states = Set() for e ∈ edges(cl_h) - if e in edge_states continue end + if e in edge_states + continue + end v, w = src(e), dst(e) v1, v2, _ = v w1, w2, _ = w - if (v1, v2) == (w1, w2) continue end + if (v1, v2) == (w1, w2) + continue + end add_edge!(new_cl_h, (v1, v2), (w1, w2)) @@ -324,8 +339,12 @@ the provided temperature parameter `β`. The merged energy values, left projector `pl`, and right projector `pr` are computed based on the interactions between the original vertices and their respective projectors. """ -function merge_vertices_cl_h(cl_h::LabelledGraph{S, T}, β::Real, node1::NTuple{3, Int64}, node2::NTuple{3, Int64} - ) where {S, T} +function merge_vertices_cl_h( + cl_h::LabelledGraph{S,T}, + β::Real, + node1::NTuple{3,Int64}, + node2::NTuple{3,Int64}, +) where {S,T} i1, j1, _ = node1 i2, j2, _ = node2 @@ -378,7 +397,7 @@ If the vertex exists in the graph and has associated energy values, it returns t The local energy values are typically obtained from the spectrum associated with the vertex. """ -function local_energy(cl_h::LabelledGraph{S, T}, v::NTuple{3, Int64}) where {S, T} +function local_energy(cl_h::LabelledGraph{S,T}, v::NTuple{3,Int64}) where {S,T} has_vertex(cl_h, v) ? get_prop(cl_h, v, :spectrum).energies : zeros(1) end @@ -401,7 +420,11 @@ if there is a directed edge from `v` to `w`, it returns the transpose of the ene otherwise, it returns a matrix of zeros. The interaction energy values represent the energy associated with the interaction or connection between the two vertices. """ -function interaction_energy(cl_h::LabelledGraph{S, T}, v::NTuple{3, Int64}, w::NTuple{3, Int64}) where {S, T} +function interaction_energy( + cl_h::LabelledGraph{S,T}, + v::NTuple{3,Int64}, + w::NTuple{3,Int64}, +) where {S,T} if has_edge(cl_h, w, v) get_prop(cl_h, w, v, :en)' elseif has_edge(cl_h, v, w) @@ -430,7 +453,11 @@ If there is a directed edge from `w` to `v`, it returns the index of right proje if there is a directed edge from `v` to `w`, it returns the index of left projector (`:ipl`). If no edge exists between the vertices, it returns a vector of ones. """ -function projector(cl_h::LabelledGraph{S, T}, v::NTuple{N, Int64}, w::NTuple{N, Int64}) where {S, T, N} +function projector( + cl_h::LabelledGraph{S,T}, + v::NTuple{N,Int64}, + w::NTuple{N,Int64}, +) where {S,T,N} if has_edge(cl_h, w, v) idx_p = get_prop(cl_h, w, v, :ipr) p = get_projector!(get_prop(cl_h, :pool_of_projectors), idx_p, :CPU) @@ -438,17 +465,20 @@ function projector(cl_h::LabelledGraph{S, T}, v::NTuple{N, Int64}, w::NTuple{N, idx_p = get_prop(cl_h, v, w, :ipl) p = get_projector!(get_prop(cl_h, :pool_of_projectors), idx_p, :CPU) else - p = ones(Int, v ∈ vertices(cl_h) ? length(get_prop(cl_h, v, :spectrum).energies) : 1) + p = ones( + Int, + v ∈ vertices(cl_h) ? length(get_prop(cl_h, v, :spectrum).energies) : 1, + ) end end -function fuse_projectors(projectors::NTuple{N, K}) where {N, K} +function fuse_projectors(projectors::NTuple{N,K}) where {N,K} fused, transitions_matrix = rank_reveal(hcat(projectors...), :PE) transitions = Tuple(Array(t) for t ∈ eachcol(transitions_matrix)) fused, transitions end -function outer_projector(p1::Array{T, 1}, p2::Array{T, 1}) where T <: Number +function outer_projector(p1::Array{T,1}, p2::Array{T,1}) where {T<:Number} reshape(reshape(p1, :, 1) .+ maximum(p1) .* reshape(p2 .- 1, 1, :), :) end @@ -468,10 +498,10 @@ This constructor function creates a sparse column-compressed (CSC) matrix of ele column indices `p` and values. The resulting matrix has non-zero values at the specified column indices, while all other elements are zero. The `SparseCSC` constructor is useful for creating sparse matrices with specific column indices and values efficiently. """ -function SparseCSC(::Type{R}, p::Vector{Int64}) where R <: Real +function SparseCSC(::Type{R}, p::Vector{Int64}) where {R<:Real} n = length(p) mp = maximum(p) cn = collect(1:n) co = ones(R, n) sparse(p, cn, co, mp, n) -end \ No newline at end of file +end diff --git a/src/clustered_hamiltonian.jl b/src/clustered_hamiltonian.jl index 996c886..f5c3f78 100644 --- a/src/clustered_hamiltonian.jl +++ b/src/clustered_hamiltonian.jl @@ -1,5 +1,4 @@ -export - clustered_hamiltonian, +export clustered_hamiltonian, rank_reveal, split_into_clusters, decode_clustered_hamiltonian_state, @@ -31,9 +30,11 @@ Each cluster is represented by a vertex from the Ising graph. The `split_into_clusters` function is useful for organizing and analyzing spins in complex spin systems, particularly in the context of clustered Hamiltonian. """ -function split_into_clusters(ig::LabelledGraph{G, L}, assignment_rule) where {G, L} +function split_into_clusters(ig::LabelledGraph{G,L}, assignment_rule) where {G,L} cluster_id_to_verts = Dict(i => L[] for i in values(assignment_rule)) - for v in vertices(ig) push!(cluster_id_to_verts[assignment_rule[v]], v) end + for v in vertices(ig) + push!(cluster_id_to_verts[assignment_rule[v]], v) + end Dict(i => first(cluster(ig, verts)) for (i, verts) ∈ cluster_id_to_verts) end @@ -61,11 +62,16 @@ and a cluster assignment rule, which maps Ising graph vertices to clusters. function clustered_hamiltonian( ig::IsingGraph, num_states_cl::Int; - spectrum::Function=full_spectrum, - cluster_assignment_rule::Dict{Int, L} # e.g. square lattice -) where L + spectrum::Function = full_spectrum, + cluster_assignment_rule::Dict{Int,L}, # e.g. square lattice +) where {L} ns = Dict(i => num_states_cl for i ∈ Set(values(cluster_assignment_rule))) - clustered_hamiltonian(ig, ns, spectrum=spectrum, cluster_assignment_rule=cluster_assignment_rule) + clustered_hamiltonian( + ig, + ns, + spectrum = spectrum, + cluster_assignment_rule = cluster_assignment_rule, + ) end """ @@ -92,18 +98,16 @@ and a cluster assignment rule, which maps Ising graph vertices to clusters. """ function clustered_hamiltonian( ig::IsingGraph, - num_states_cl::Dict{T, Int}; - spectrum::Function=full_spectrum, - cluster_assignment_rule::Dict{Int, T} -) where T - cl_h = LabelledGraph{MetaDiGraph}( - sort(unique(values(cluster_assignment_rule))) - ) + num_states_cl::Dict{T,Int}; + spectrum::Function = full_spectrum, + cluster_assignment_rule::Dict{Int,T}, +) where {T} + cl_h = LabelledGraph{MetaDiGraph}(sort(unique(values(cluster_assignment_rule)))) lp = PoolOfProjectors{Int}() for (v, cl) ∈ split_into_clusters(ig, cluster_assignment_rule) - sp = spectrum(cl, num_states=get(num_states_cl, v, basis_size(cl))) + sp = spectrum(cl, num_states = get(num_states_cl, v, basis_size(cl))) set_props!(cl_h, v, Dict(:cluster => cl, :spectrum => sp)) end @@ -112,8 +116,8 @@ function clustered_hamiltonian( outer_edges, J = inter_cluster_edges(ig, cl1, cl2) if !isempty(outer_edges) - ind1 = any(i -> i != 0, J, dims=2) - ind2 = any(i -> i != 0, J, dims=1) + ind1 = any(i -> i != 0, J, dims = 2) + ind2 = any(i -> i != 0, J, dims = 1) ind1 = reshape(ind1, length(ind1)) ind2 = reshape(ind2, length(ind2)) JJ = J[ind1, ind2] @@ -129,7 +133,10 @@ function clustered_hamiltonian( add_edge!(cl_h, v, w) set_props!( - cl_h, v, w, Dict(:outer_edges => outer_edges, :ipl => ipl, :en => en, :ipr => ipr) + cl_h, + v, + w, + Dict(:outer_edges => outer_edges, :ipl => ipl, :en => en, :ipr => ipr), ) end end @@ -160,11 +167,16 @@ If you want to specify custom cluster sizes, use the alternative version of this passing a `Dict{T, Int}` containing the number of states per cluster as `num_states_cl`. """ function clustered_hamiltonian( - ig::IsingGraph; - spectrum::Function=full_spectrum, - cluster_assignment_rule::Dict{Int, T} - ) where T - clustered_hamiltonian(ig, Dict{T, Int}(), spectrum=spectrum, cluster_assignment_rule=cluster_assignment_rule) + ig::IsingGraph; + spectrum::Function = full_spectrum, + cluster_assignment_rule::Dict{Int,T}, +) where {T} + clustered_hamiltonian( + ig, + Dict{T,Int}(), + spectrum = spectrum, + cluster_assignment_rule = cluster_assignment_rule, + ) end # """ @@ -216,8 +228,11 @@ returns a dictionary mapping each Ising graph vertex to its corresponding spin v This function assumes that the state has the same order as the vertices in the clustered Hamiltonian. It decodes the state consistently based on the cluster assignments and spectra of the clustered Hamiltonian. """ -function decode_clustered_hamiltonian_state(cl_h::LabelledGraph{S, T}, state::Vector{Int}) where {S, T} - ret = Dict{Int, Int}() +function decode_clustered_hamiltonian_state( + cl_h::LabelledGraph{S,T}, + state::Vector{Int}, +) where {S,T} + ret = Dict{Int,Int}() for (i, vert) ∈ zip(state, vertices(cl_h)) spins = get_prop(cl_h, vert, :cluster).labels states = get_prop(cl_h, vert, :spectrum).states @@ -248,9 +263,11 @@ This function computes the energy by summing the energies associated with indivi clusters and the interaction energies between clusters. It takes into account the cluster spectra and projectors stored in the clustered Hamiltonian. """ -function energy(cl_h::LabelledGraph{S, T}, σ::Dict{T, Int}) where {S, T} +function energy(cl_h::LabelledGraph{S,T}, σ::Dict{T,Int}) where {S,T} en_cl_h = 0.0 - for v ∈ vertices(cl_h) en_cl_h += get_prop(cl_h, v, :spectrum).energies[σ[v]] end + for v ∈ vertices(cl_h) + en_cl_h += get_prop(cl_h, v, :spectrum).energies[σ[v]] + end for edge ∈ edges(cl_h) idx_pl = get_prop(cl_h, edge, :ipl) pl = get_projector!(get_prop(cl_h, :pool_of_projectors), idx_pl, :CPU) @@ -281,7 +298,7 @@ The function checks if there is an interaction edge between the two sites (i, j) If such edges exist, it retrieves the interaction energy matrix, projectors, and calculates the interaction energy. If no interaction edge is found, it returns a zero matrix. """ -function energy_2site(cl_h::LabelledGraph{S, T}, i::Int, j::Int) where {S, T} +function energy_2site(cl_h::LabelledGraph{S,T}, i::Int, j::Int) where {S,T} # matrix of interaction energies between two nodes if has_edge(cl_h, (i, j, 1), (i, j, 2)) en12 = copy(get_prop(cl_h, (i, j, 1), (i, j, 2), :en)) @@ -324,22 +341,18 @@ If such edges exist, it retrieves the bond energy matrix and projectors and calc If no bond edge is found, it returns a zero vector. """ function bond_energy( - cl_h::LabelledGraph{S, T}, - cl_h_u::NTuple{N, Int64}, - cl_h_v::NTuple{N, Int64}, - σ::Int - ) where {S, T, N} + cl_h::LabelledGraph{S,T}, + cl_h_u::NTuple{N,Int64}, + cl_h_v::NTuple{N,Int64}, + σ::Int, +) where {S,T,N} if has_edge(cl_h, cl_h_u, cl_h_v) - ipu, en, ipv = get_prop.( - Ref(cl_h), Ref(cl_h_u), Ref(cl_h_v), (:ipl, :en, :ipr) - ) + ipu, en, ipv = get_prop.(Ref(cl_h), Ref(cl_h_u), Ref(cl_h_v), (:ipl, :en, :ipr)) pu = get_projector!(get_prop(cl_h, :pool_of_projectors), ipu, :CPU) pv = get_projector!(get_prop(cl_h, :pool_of_projectors), ipv, :CPU) @inbounds energies = en[pu, pv[σ]] elseif has_edge(cl_h, cl_h_v, cl_h_u) - ipv, en, ipu = get_prop.( - Ref(cl_h), Ref(cl_h_v), Ref(cl_h_u), (:ipl, :en, :ipr) - ) + ipv, en, ipu = get_prop.(Ref(cl_h), Ref(cl_h_v), Ref(cl_h_u), (:ipl, :en, :ipr)) pu = get_projector!(get_prop(cl_h, :pool_of_projectors), ipu, :CPU) pv = get_projector!(get_prop(cl_h, :pool_of_projectors), ipv, :CPU) @inbounds energies = en[pv[σ], pu] @@ -364,7 +377,7 @@ This function returns the size (number of states) of a cluster in a clustered Ha The function retrieves the spectrum associated with the specified cluster and returns the length of the energy vector in that spectrum. """ -function cluster_size(clustered_hamiltonian::LabelledGraph{S, T}, vertex::T) where {S, T} +function cluster_size(clustered_hamiltonian::LabelledGraph{S,T}, vertex::T) where {S,T} length(get_prop(clustered_hamiltonian, vertex, :spectrum).energies) end @@ -387,7 +400,11 @@ The function generates all possible states for the clusters in the clustered Ham calculates their energies, and computes the probability distribution based on the given inverse temperature parameter. It then calculates the conditional probability of the specified target state by summing the probabilities of states that match the target state. """ -function exact_cond_prob(clustered_hamiltonian::LabelledGraph{S, T}, beta, target_state::Dict) where {S, T} +function exact_cond_prob( + clustered_hamiltonian::LabelledGraph{S,T}, + beta, + target_state::Dict, +) where {S,T} # TODO: Not going to work without PoolOfProjectors ver = vertices(clustered_hamiltonian) rank = cluster_size.(Ref(clustered_hamiltonian), ver) @@ -418,7 +435,7 @@ It then updates the spectrum of each cluster in `new_cl_h` by selecting the spec Additionally, it updates the interactions and projectors between clusters based on the retained states. The resulting `new_cl_h` represents a truncated version of the original Hamiltonian. """ -function truncate_clustered_hamiltonian(cl_h::LabelledGraph{S, T}, states::Dict) where {S, T} +function truncate_clustered_hamiltonian(cl_h::LabelledGraph{S,T}, states::Dict) where {S,T} new_cl_h = LabelledGraph{MetaDiGraph}(vertices(cl_h)) new_lp = PoolOfProjectors{Int}() @@ -427,7 +444,7 @@ function truncate_clustered_hamiltonian(cl_h::LabelledGraph{S, T}, states::Dict) cl = get_prop(cl_h, v, :cluster) sp = get_prop(cl_h, v, :spectrum) if sp.states == Vector{Int64}[] - sp = Spectrum(sp.energies[states[v]], sp.states, [1,]) + sp = Spectrum(sp.energies[states[v]], sp.states, [1]) else sp = Spectrum(sp.energies[states[v]], sp.states[states[v]]) end @@ -451,17 +468,22 @@ function truncate_clustered_hamiltonian(cl_h::LabelledGraph{S, T}, states::Dict) ipl = add_projector!(new_lp, pl_transition) ipr = add_projector!(new_lp, pr_transition) set_props!( - new_cl_h, v, w, Dict(:outer_edges => outer_edges, :ipl => ipl, :en => en, :ipr => ipr) - ) + new_cl_h, + v, + w, + Dict(:outer_edges => outer_edges, :ipl => ipl, :en => en, :ipr => ipr), + ) end set_props!(new_cl_h, Dict(:pool_of_projectors => new_lp)) new_cl_h end -function clustered_hamiltonian(fname::String, - Nx::Union{Integer, Nothing}=nothing, - Ny::Union{Integer, Nothing}=nothing) +function clustered_hamiltonian( + fname::String, + Nx::Union{Integer,Nothing} = nothing, + Ny::Union{Integer,Nothing} = nothing, +) loaded_rmf = load_openGM(fname, Nx, Ny) functions = loaded_rmf["fun"] @@ -481,20 +503,33 @@ function clustered_hamiltonian(fname::String, y, x = index Eng = functions[value]' sp = Spectrum(collect(Eng), Vector{Vector{Int}}[], zeros(Int, N[y+1, x+1])) - set_props!(cl_h, (x+1, y+1), Dict(:spectrum => sp)) + set_props!(cl_h, (x + 1, y + 1), Dict(:spectrum => sp)) elseif length(index) == 4 y1, x1, y2, x2 = index add_edge!(cl_h, (x1 + 1, y1 + 1), (x2 + 1, y2 + 1)) Eng = functions[value] ipl = add_projector!(lp, collect(1:N[y1+1, x1+1])) ipr = add_projector!(lp, collect(1:N[y2+1, x2+1])) - set_props!(cl_h, (x1 + 1, y1 + 1), (x2 + 1, y2 + 1), Dict(:outer_edges=> ((x1 + 1, y1 + 1), (x2 + 1, y2 + 1)), - :en => Eng, :ipl => ipl, :ipr => ipr)) + set_props!( + cl_h, + (x1 + 1, y1 + 1), + (x2 + 1, y2 + 1), + Dict( + :outer_edges => ((x1 + 1, y1 + 1), (x2 + 1, y2 + 1)), + :en => Eng, + :ipl => ipl, + :ipr => ipr, + ), + ) else - throw(ErrorException("Something is wrong with factor index, it has length $(length(index))")) + throw( + ErrorException( + "Something is wrong with factor index, it has length $(length(index))", + ), + ) end end - + set_props!(cl_h, Dict(:pool_of_projectors => lp, :Nx => X, :Ny => Y)) cl_h -end \ No newline at end of file +end diff --git a/src/ising.jl b/src/ising.jl index 4229ec3..fcdb63e 100644 --- a/src/ising.jl +++ b/src/ising.jl @@ -1,7 +1,6 @@ using LabelledGraphs -export - ising_graph, +export ising_graph, rank_vec, cluster, rank, @@ -13,8 +12,8 @@ export prune, inter_cluster_edges -const Instance = Union{String, Dict} -const IsingGraph{T} = LabelledGraph{MetaGraph{Int, T}} +const Instance = Union{String,Dict} +const IsingGraph{T} = LabelledGraph{MetaGraph{Int,T}} function unique_nodes(ising_tuples) sort(collect(Set(Iterators.flatten((i, j) for (i, j, _) ∈ ising_tuples)))) @@ -43,9 +42,14 @@ It assigns interaction strengths to edges between spins and optionally scales th The `rank_override` dictionary can be used to specify the rank (number of states) for individual vertices, allowing customization of the Ising model. Convention: H = scale * sum_{i, j} (J_{ij} * s_i * s_j + J_{ii} * s_i) """ -function ising_graph(::Type{T}, inst::Instance; scale::Real=1, rank_override::Dict=Dict{Int, Int}()) where T +function ising_graph( + ::Type{T}, + inst::Instance; + scale::Real = 1, + rank_override::Dict = Dict{Int,Int}(), +) where {T} if inst isa String - ising = CSV.File(inst, types = [Int, Int, T], header=0, comment = "#") + ising = CSV.File(inst, types = [Int, Int, T], header = 0, comment = "#") else ising = [(i, j, T(J)) for ((i, j), J) ∈ inst] end @@ -67,10 +71,10 @@ function ising_graph(::Type{T}, inst::Instance; scale::Real=1, rank_override::Di ig end -function ising_graph(inst::Instance; scale::Real=1, rank_override::Dict=Dict{Int, Int}()) +function ising_graph(inst::Instance; scale::Real = 1, rank_override::Dict = Dict{Int,Int}()) ising_graph(Float64, inst; scale = scale, rank_override = rank_override) end -Base.eltype(ig::IsingGraph{T}) where T = T +Base.eltype(ig::IsingGraph{T}) where {T} = T rank_vec(ig::IsingGraph) = Int[get_prop((ig), v, :rank) for v ∈ vertices(ig)] basis_size(ig::IsingGraph) = prod(rank_vec(ig)) @@ -92,7 +96,7 @@ The coupling strengths are represented as a matrix, where each element `(i, j)` The function iterates over the edges of the Ising graph and extracts the interaction strengths associated with each edge, populating the `J` matrix accordingly. """ -function couplings(ig::IsingGraph{T}) where T +function couplings(ig::IsingGraph{T}) where {T} J = zeros(T, nv(ig), nv(ig)) for edge ∈ edges(ig) i = ig.reverse_label_map[src(edge)] @@ -125,8 +129,13 @@ where each element `(i, j)` corresponds to the interaction strength between clus The function first identifies the outer edges that connect vertices between the two clusters in the context of the larger Ising graph `ig`. It then computes the interaction strengths associated with these outer edges and populates the dense adjacency matrix `J` accordingly. """ -function inter_cluster_edges(ig::IsingGraph{T}, cl1::IsingGraph{T}, cl2::IsingGraph{T}) where T - outer_edges = [LabelledEdge(i, j) for i ∈ vertices(cl1), j ∈ vertices(cl2) if has_edge(ig, i, j)] +function inter_cluster_edges( + ig::IsingGraph{T}, + cl1::IsingGraph{T}, + cl2::IsingGraph{T}, +) where {T} + outer_edges = + [LabelledEdge(i, j) for i ∈ vertices(cl1), j ∈ vertices(cl2) if has_edge(ig, i, j)] J = zeros(T, nv(cl1), nv(cl2)) for e ∈ outer_edges i, j = cl1.reverse_label_map[src(e)], cl2.reverse_label_map[dst(e)] @@ -154,13 +163,17 @@ magnetic field (`h`) that is not approximately equal to zero within the specifie The function returns a pruned version of the input Ising graph, where non-existing spins and their associated properties are removed. """ -function prune(ig::IsingGraph; atol::Real=1e-14) +function prune(ig::IsingGraph; atol::Real = 1e-14) to_keep = vcat( findall(!iszero, degree(ig)), - findall(x -> iszero(degree(ig, x)) && !isapprox(get_prop(ig, x, :h), 0, atol=atol), vertices(ig)) + findall( + x -> + iszero(degree(ig, x)) && !isapprox(get_prop(ig, x, :h), 0, atol = atol), + vertices(ig), + ), ) gg = ig[ig.labels[to_keep]] labels = collect(vertices(gg.inner_graph)) - reverse_label_map = Dict(i => i for i=1:nv(gg.inner_graph)) + reverse_label_map = Dict(i => i for i = 1:nv(gg.inner_graph)) LabelledGraph(labels, gg.inner_graph, reverse_label_map) end diff --git a/src/lattice.jl b/src/lattice.jl index 5ef94c4..1b90ff3 100644 --- a/src/lattice.jl +++ b/src/lattice.jl @@ -1,5 +1,4 @@ -export - super_square_lattice, +export super_square_lattice, pegasus_lattice, pegasus_lattice_masoud, pegasus_lattice_tomek, @@ -28,10 +27,10 @@ where m is the number of columns, n is the number of rows and t denotes the numb The `size` tuple represents the dimensions of the super square lattice. The function creates a dictionary where ising graph coordinates are associated with their corresponding lattice coordinates. """ -function super_square_lattice(size::NTuple{5, Int}) +function super_square_lattice(size::NTuple{5,Int}) m, um, n, un, t = size old = LinearIndices((1:t, 1:un, 1:n, 1:um, 1:m)) - Dict(old[k, uj, j, ui, i] => (i, j) for i=1:m, ui=1:um, j=1:n, uj=1:un, k=1:t) + Dict(old[k, uj, j, ui, i] => (i, j) for i = 1:m, ui = 1:um, j = 1:n, uj = 1:un, k = 1:t) end """ @@ -52,12 +51,12 @@ n is the number of rows and t denotes the number of spins stored in the cluster. The `size` tuple represents the dimensions of the simplified super square lattice. The function internally adds the required dimensions `(1, 1)` to make it compatible with the `super_square_lattice` function, which deals with five dimensions. """ -function super_square_lattice(size::NTuple{3, Int}) +function super_square_lattice(size::NTuple{3,Int}) m, n, t = size super_square_lattice((m, 1, n, 1, t)) end -pegasus_lattice(size::NTuple{2, Int}) = pegasus_lattice((size[1], size[2], 3)) +pegasus_lattice(size::NTuple{2,Int}) = pegasus_lattice((size[1], size[2], 3)) """ $(TYPEDSIGNATURES) @@ -75,64 +74,82 @@ based on the specified size of the Pegasus lattice in three dimensions: `(m, n, The `pegasus_lattice` allows you to build the graph relevant for D-Wave Pegasus architecture. """ -function pegasus_lattice(size::NTuple{3, Int}) +function pegasus_lattice(size::NTuple{3,Int}) m, n, t = size old = LinearIndices((1:8*t, 1:n, 1:m)) map = Dict( - old[k, j, i] => (i, j, 1) for i=1:m, j=1:n, k ∈ (p * 8 + q for p ∈ 0 : t-1, q ∈ 1:4) + old[k, j, i] => (i, j, 1) for i = 1:m, j = 1:n, + k ∈ (p * 8 + q for p ∈ 0:t-1, q ∈ 1:4) ) - for i=1:m, j=1:n, k ∈ (p * 8 + q for p ∈ 0 : t-1, q ∈ 5:8) + for i = 1:m, j = 1:n, k ∈ (p * 8 + q for p ∈ 0:t-1, q ∈ 5:8) push!(map, old[k, j, i] => (i, j, 2)) end map end # TODO masoud / tomek should be removed from function names -function pegasus_lattice_alternative(size::NTuple{3, Int}) +function pegasus_lattice_alternative(size::NTuple{3,Int}) m, n, t = size old = LinearIndices((1:8*t, 1:n, 1:m)) map = Dict( - old[k, j, i] => (i, j, 2) for i=1:m, j=1:n, k ∈ (p * 8 + q for p ∈ 0 : t-1, q ∈ 1:4) + old[k, j, i] => (i, j, 2) for i = 1:m, j = 1:n, + k ∈ (p * 8 + q for p ∈ 0:t-1, q ∈ 1:4) ) - for i=1:m, j=1:n, k ∈ (p * 8 + q for p ∈ 0 : t-1, q ∈ 5:8) + for i = 1:m, j = 1:n, k ∈ (p * 8 + q for p ∈ 0:t-1, q ∈ 5:8) push!(map, old[k, j, i] => (i, j, 1)) end map end -function pegasus_lattice_old_numering(size::NTuple{3, Int}) +function pegasus_lattice_old_numering(size::NTuple{3,Int}) m, n, t = size old = LinearIndices((1:8*t, 1:n, 1:m)) map = Dict( - old[k, j, i] => (i, n-j+1, 2) for i=1:m, j=1:n, k ∈ (p * 8 + q for p ∈ 0 : t-1, q ∈ 1:4) + old[k, j, i] => (i, n - j + 1, 2) for i = 1:m, j = 1:n, + k ∈ (p * 8 + q for p ∈ 0:t-1, q ∈ 1:4) ) - for i=1:m, j=1:n, k ∈ (p * 8 + q for p ∈ 0 : t-1, q ∈ 5:8) - push!(map, old[k, j, i] => (i, n-j+1, 1)) + for i = 1:m, j = 1:n, k ∈ (p * 8 + q for p ∈ 0:t-1, q ∈ 5:8) + push!(map, old[k, j, i] => (i, n - j + 1, 1)) end map end -function zephyr_lattice_z1(size::NTuple{3, Int}) +function zephyr_lattice_z1(size::NTuple{3,Int}) m, n, t = size # t is identical to dwave (Tile parameter for the Zephyr lattice) - map = Dict{Int, NTuple{3, Int}}() - for i=1:2*n, j ∈ 1:2*m + map = Dict{Int,NTuple{3,Int}}() + for i = 1:2*n, j ∈ 1:2*m for p in p_func(i, j, t, n, m) - push!(map, (i-1)*(2*n*t) + (j-1)*(2*m*t) + p*n + (i-1)*(j%2) + 1 => (i, j, 1)) + push!( + map, + (i - 1) * (2 * n * t) + + (j - 1) * (2 * m * t) + + p * n + + (i - 1) * (j % 2) + + 1 => (i, j, 1), + ) end for q ∈ q_func(i, j, t, n, m) - push!(map, 2*t*(2*n+1) + (i-1)*(2*n*t) + (j%2)*(2*m*t) + q*m + (j-1)*(i-1) + 1 => (i, j, 2)) + push!( + map, + 2 * t * (2 * n + 1) + + (i - 1) * (2 * n * t) + + (j % 2) * (2 * m * t) + + q * m + + (j - 1) * (i - 1) + + 1 => (i, j, 2), + ) end end map end function j_function(i::Int, n::Int) - i ∈ collect(1:n) && return collect((n + 1 - i):(n + i)) - collect((i-n):(3*n + 1 - i)) + i ∈ collect(1:n) && return collect((n+1-i):(n+i)) + collect((i-n):(3*n+1-i)) end -zephyr_lattice(size::NTuple{2, Int}) = zephyr_lattice((size[1], size[2], 4)) +zephyr_lattice(size::NTuple{2,Int}) = zephyr_lattice((size[1], size[2], 4)) """ $(TYPEDSIGNATURES) @@ -150,48 +167,58 @@ coordinates based on the specified size of the Zephyr lattice in three dimension The `zephyr_lattice` allows you to build the graph relevant for D-Wave Zephyr architecture. """ -function zephyr_lattice(size::NTuple{3, Int}) +function zephyr_lattice(size::NTuple{3,Int}) m, n, t = size - zephyr_lattice_5tuple_rotated(m+1, n+1, zephyr_lattice_5tuple((Int(m/2), Int(n/2), t))) + zephyr_lattice_5tuple_rotated( + m + 1, + n + 1, + zephyr_lattice_5tuple((Int(m / 2), Int(n / 2), t)), + ) end -function zephyr_lattice_5tuple(size::NTuple{3, Int}) - m, n , t = size # t is identical to dwave (Tile parameter for the Zephyr lattice) - map = Dict{Int, NTuple{3, Int}}() +function zephyr_lattice_5tuple(size::NTuple{3,Int}) + m, n, t = size # t is identical to dwave (Tile parameter for the Zephyr lattice) + map = Dict{Int,NTuple{3,Int}}() # ( u, w, k, ζ, z) - for u = 0 - for w ∈ 0:2:2*m, k ∈ 0:t-1, ζ ∈ 0:1, (i,z) ∈ enumerate(0:n-1) - push!(map, zephyr_to_linear(m, t, (u,w,k,ζ,z)) => (2*i, w + 1, 1)) + for u in 0 + for w ∈ 0:2:2*m, k ∈ 0:t-1, ζ ∈ 0:1, (i, z) ∈ enumerate(0:n-1) + push!(map, zephyr_to_linear(m, t, (u, w, k, ζ, z)) => (2 * i, w + 1, 1)) end for w ∈ 1:2:2*m, k ∈ 0:t-1, ζ ∈ 0:1, z ∈ 0:n-1 - push!(map, zephyr_to_linear(m, t, (u,w,k,ζ,z)) => (2*z + 2*ζ + 1, w + 1, 1)) + push!( + map, + zephyr_to_linear(m, t, (u, w, k, ζ, z)) => (2 * z + 2 * ζ + 1, w + 1, 1), + ) end end - for u = 1 - for w ∈ 0:2:2*m, k ∈ 0:t-1, ζ ∈ 0:1, (i,z) ∈ enumerate(0:n-1) - push!(map, zephyr_to_linear(m, t, (u,w,k,ζ,z)) => (w + 1, 2*i, 2)) + for u in 1 + for w ∈ 0:2:2*m, k ∈ 0:t-1, ζ ∈ 0:1, (i, z) ∈ enumerate(0:n-1) + push!(map, zephyr_to_linear(m, t, (u, w, k, ζ, z)) => (w + 1, 2 * i, 2)) end for w ∈ 1:2:2*m, k ∈ 0:t-1, ζ ∈ 0:1, z ∈ 0:n-1 - push!(map, zephyr_to_linear(m, t, (u,w,k,ζ,z)) => (w + 1, 2*z + 2*ζ + u, 2)) + push!( + map, + zephyr_to_linear(m, t, (u, w, k, ζ, z)) => (w + 1, 2 * z + 2 * ζ + u, 2), + ) end end map end function rotate(m::Int, n::Int) - new_dict = Dict{NTuple{3, Int}, NTuple{3, Int}}() + new_dict = Dict{NTuple{3,Int},NTuple{3,Int}}() for (k, j) ∈ enumerate(1:2:m) - for (l,i) ∈ enumerate(n-1:-2:1) - push!(new_dict, (i,j,1) => (i/2 + k - 1, l + k - 1, 1)) - push!(new_dict, (i,j,2) => (i/2 + k - 1, l + k - 1, 2)) + for (l, i) ∈ enumerate(n-1:-2:1) + push!(new_dict, (i, j, 1) => (i / 2 + k - 1, l + k - 1, 1)) + push!(new_dict, (i, j, 2) => (i / 2 + k - 1, l + k - 1, 2)) end end for (k, j) ∈ enumerate(2:2:m) - for (l,i) ∈ enumerate(n:-2:1) - push!(new_dict, (i,j,1) => ((i-1)/2 + k, l + k - 1, 1)) - push!(new_dict, (i,j,2) => ((i-1)/2 + k, l + k - 1, 2)) + for (l, i) ∈ enumerate(n:-2:1) + push!(new_dict, (i, j, 1) => ((i - 1) / 2 + k, l + k - 1, 1)) + push!(new_dict, (i, j, 2) => ((i - 1) / 2 + k, l + k - 1, 2)) end end new_dict @@ -207,16 +234,16 @@ function empty_clusters(m::Int, n::Int) (count, reverse(ii)) end -function zephyr_lattice_5tuple_rotated(m::Int, n::Int, map::Dict{Int, NTuple{3, Int}}) +function zephyr_lattice_5tuple_rotated(m::Int, n::Int, map::Dict{Int,NTuple{3,Int}}) rotated_map = rotate(m, n) - new_map = Dict{Int, NTuple{3, Int}}() + new_map = Dict{Int,NTuple{3,Int}}() (empty, ii) = empty_clusters(m, n) for k in keys(map) push!(new_map, k => rotated_map[map[k]]) end - empty_vertices = empty_indexing(m,n) - for (k,l) ∈ enumerate(empty_vertices) + empty_vertices = empty_indexing(m, n) + for (k, l) ∈ enumerate(empty_vertices) push!(new_map, -k => l) end new_map @@ -224,7 +251,7 @@ end function empty_indexing(m::Int, n::Int) (empty, ii) = empty_clusters(m, n) - p = Int((m-1)/2) + p = Int((m - 1) / 2) empty_vertices = [] for (k, l) ∈ enumerate(ii) for i ∈ 1:l @@ -234,7 +261,7 @@ function empty_indexing(m::Int, n::Int) push!(empty_vertices, (k, i + m - p + k - 1, 2)) end end - for (k,l) ∈ enumerate(reverse(ii)) + for (k, l) ∈ enumerate(reverse(ii)) for i ∈ 1:l push!(empty_vertices, (k + m - p, i, 1)) push!(empty_vertices, (k + m - p, i, 2)) @@ -245,24 +272,24 @@ function empty_indexing(m::Int, n::Int) empty_vertices end -function periodic_lattice(size::NTuple{3, Int}) +function periodic_lattice(size::NTuple{3,Int}) mm, nn, tt = size - m, n = 2*mm, 2*nn + m, n = 2 * mm, 2 * nn map = super_square_lattice((m, n, 1)) - new_map = Dict{Int, NTuple{2, Int}}() + new_map = Dict{Int,NTuple{2,Int}}() for (key, val) ∈ map i, j = val - if i <= m/2 - if j <= m/2 + if i <= m / 2 + if j <= m / 2 push!(new_map, key => (i, j)) - elseif j > m/2 - push!(new_map, key => (i, m-j+1)) + elseif j > m / 2 + push!(new_map, key => (i, m - j + 1)) end - elseif i > m/2 - if j <= m/2 - push!(new_map, key => (m-i+1, j)) - elseif j > m/2 - push!(new_map, key => (m-i+1, m-j+1)) + elseif i > m / 2 + if j <= m / 2 + push!(new_map, key => (m - i + 1, j)) + elseif j > m / 2 + push!(new_map, key => (m - i + 1, m - j + 1)) end end end diff --git a/src/spectrum.jl b/src/spectrum.jl index d38df75..8efabfa 100644 --- a/src/spectrum.jl +++ b/src/spectrum.jl @@ -1,5 +1,4 @@ -export - all_states, +export all_states, local_basis, gibbs_tensor, brute_force, @@ -12,7 +11,7 @@ export @inline idx(σ::Int) = (σ == -1) ? 1 : σ + 1 @inline local_basis(d::Int) = union(-1, 1:d-1) -all_states(rank::Union{Vector, NTuple}) = Iterators.product(local_basis.(rank)...) +all_states(rank::Union{Vector,NTuple}) = Iterators.product(local_basis.(rank)...) const State = Vector{Int} @@ -58,10 +57,10 @@ This function takes a matrix of binary vectors, where each row represents a bina # Returns: - `Vector{Int}`: An array of integer representations of the binary vectors. """ -function matrix_to_integers(matrix::Vector{Vector{T}}) where T +function matrix_to_integers(matrix::Vector{Vector{T}}) where {T} nrows = length(matrix[1]) multipliers = 2 .^ collect(0:nrows-1) - div.((hcat(matrix...)' .+ 1) , 2) * multipliers + div.((hcat(matrix...)' .+ 1), 2) * multipliers end """ @@ -99,15 +98,15 @@ The energy is calculated based on the interactions between spins and their assoc # Returns: - `T`: The energy of the state in the Ising graph. """ -function energy(ig::IsingGraph{T}, ig_state::Dict{Int, Int}) where T +function energy(ig::IsingGraph{T}, ig_state::Dict{Int,Int}) where {T} en = zero(T) for (i, σ) ∈ ig_state en += get_prop(ig, i, :h) * σ for (j, η) ∈ ig_state if has_edge(ig, i, j) - en += T(1/2) * σ * get_prop(ig, i, j, :J) * η + en += T(1 / 2) * σ * get_prop(ig, i, j, :J) * η elseif has_edge(ig, j, i) - en += T(1/2) * σ * get_prop(ig, j, i, :J) * η + en += T(1 / 2) * σ * get_prop(ig, j, i, :J) * η end end end @@ -128,7 +127,7 @@ The energy spectrum represents all possible energy levels and their associated s # Returns: - `Spectrum`: An instance of the `Spectrum` type containing the energy levels and states. """ -function Spectrum(ig::IsingGraph{T}) where T +function Spectrum(ig::IsingGraph{T}) where {T} L = nv(ig) N = 2^L energies = zeros(T, N) @@ -136,7 +135,7 @@ function Spectrum(ig::IsingGraph{T}) where T J, h = couplings(ig), biases(ig) Threads.@threads for i = 0:N-1 - σ = 2 .* digits(i, base=2, pad=L) .- 1 + σ = 2 .* digits(i, base = 2, pad = L) .- 1 @inbounds energies[i+1] = dot(σ, J, σ) + dot(h, σ) @inbounds states[i+1] = σ end @@ -158,13 +157,13 @@ The Gibbs tensor represents the conditional probabilities of states given the in # Returns: - `Matrix{T}`: A matrix representing the Gibbs tensor with conditional probabilities. """ -function gibbs_tensor(ig::IsingGraph{T}, β::T=1) where T +function gibbs_tensor(ig::IsingGraph{T}, β::T = 1) where {T} σ = collect.(all_states(rank_vec(ig))) ρ = exp.(-β .* energy(σ, ig)) ρ ./ sum(ρ) end -function brute_force(ig::IsingGraph, s::Symbol=:CPU; num_states::Int=1) +function brute_force(ig::IsingGraph, s::Symbol = :CPU; num_states::Int = 1) brute_force(ig, Val(s); num_states) end @@ -186,7 +185,7 @@ The calculation is done using brute-force enumeration, making it feasible only f # Returns: - `Spectrum`: A `Spectrum` object containing the lowest-energy states and their energies. """ -function brute_force(ig::IsingGraph{T}, ::Val{:CPU}; num_states::Int=1) where T +function brute_force(ig::IsingGraph{T}, ::Val{:CPU}; num_states::Int = 1) where {T} L = nv(ig) L == 0 && return Spectrum(zeros(T, 1), Vector{Vector{Int}}[], zeros(T, 1)) sp = Spectrum(ig) @@ -195,7 +194,7 @@ function brute_force(ig::IsingGraph{T}, ::Val{:CPU}; num_states::Int=1) where T Spectrum(sp.energies[idx], sp.states[idx]) end -function full_spectrum(ig::IsingGraph{T}; num_states::Int=1) where T +function full_spectrum(ig::IsingGraph{T}; num_states::Int = 1) where {T} nv(ig) == 0 && return Spectrum(zeros(T, 1), Vector{Vector{Int}}[], zeros(T, 1)) ig_rank = rank_vec(ig) num_states = min(num_states, prod(ig_rank)) @@ -204,6 +203,10 @@ function full_spectrum(ig::IsingGraph{T}; num_states::Int=1) where T Spectrum(energies[begin:num_states], σ[begin:num_states]) end -function inter_cluster_energy(cl1_states::Vector{State}, J::Matrix{<:Real}, cl2_states::Vector{State}) +function inter_cluster_energy( + cl1_states::Vector{State}, + J::Matrix{<:Real}, + cl2_states::Vector{State}, +) hcat(collect.(cl1_states)...)' * J * hcat(collect.(cl2_states)...) end diff --git a/src/truncate.jl b/src/truncate.jl index 4883248..af6aaca 100644 --- a/src/truncate.jl +++ b/src/truncate.jl @@ -1,5 +1,4 @@ -export - truncate_clustered_hamiltonian_2site_energy, +export truncate_clustered_hamiltonian_2site_energy, truncate_clustered_hamiltonian_1site_BP, truncate_clustered_hamiltonian_2site_BP, select_numstate_best @@ -23,14 +22,14 @@ associated with a single-site cluster. It then truncates the clustered Hamiltoni - `LabelledGraph{S, T}`: A truncated clustered Hamiltonian. """ function truncate_clustered_hamiltonian_1site_BP( - cl_h::LabelledGraph{S, T}, - num_states::Int; - beta=1.0, - tol=1e-10, - iter=1 - ) where {S, T} + cl_h::LabelledGraph{S,T}, + num_states::Int; + beta = 1.0, + tol = 1e-10, + iter = 1, +) where {S,T} states = Dict() - beliefs = belief_propagation(cl_h, beta; tol=tol, iter=iter) + beliefs = belief_propagation(cl_h, beta; tol = tol, iter = iter) for node in vertices(cl_h) indices = partialsortperm(beliefs[node], 1:min(num_states, length(beliefs[node]))) push!(states, node => indices) @@ -53,11 +52,16 @@ to keep. It computes the energies for all 2-site combinations and selects the st # Returns: - `LabelledGraph{S, T}`: A truncated clustered Hamiltonian. """ -function truncate_clustered_hamiltonian_2site_energy(cl_h::LabelledGraph{S, T}, num_states::Int) where {S, T} +function truncate_clustered_hamiltonian_2site_energy( + cl_h::LabelledGraph{S,T}, + num_states::Int, +) where {S,T} # TODO: name to be clean to make it consistent with square2 and squarestar2 states = Dict() for node in vertices(cl_h) - if node in keys(states) continue end + if node in keys(states) + continue + end i, j, _ = node E1 = copy(get_prop(cl_h, (i, j, 1), :spectrum).energies) E2 = copy(get_prop(cl_h, (i, j, 2), :spectrum).energies) @@ -73,7 +77,7 @@ end function load_file(filename) if isfile(filename) - try + try load_object(string(filename)) catch e return nothing @@ -101,20 +105,24 @@ to keep. It computes the beliefs for all 2-site combinations and selects the sta - `LabelledGraph{S, T}`: A truncated clustered Hamiltonian. """ function truncate_clustered_hamiltonian_2site_BP( - cl_h::LabelledGraph{S, T}, - beliefs::Dict, + cl_h::LabelledGraph{S,T}, + beliefs::Dict, num_states::Int, result_folder::String = "results_folder", - inst::String = "inst"; - beta=1.0 - ) where {S, T} + inst::String = "inst"; + beta = 1.0, +) where {S,T} states = Dict() saved_states = load_file(joinpath(result_folder, "$(inst).jld2")) for node in vertices(cl_h) - if node in keys(states) continue end + if node in keys(states) + continue + end i, j, _ = node - sx = has_vertex(cl_h, (i, j, 1)) ? length(get_prop(cl_h, (i, j, 1), :spectrum).energies) : 1 + sx = + has_vertex(cl_h, (i, j, 1)) ? + length(get_prop(cl_h, (i, j, 1), :spectrum).energies) : 1 E = beliefs[(i, j)] ind1, ind2 = select_numstate_best(E, sx, num_states) push!(states, (i, j, 1) => ind1) @@ -163,16 +171,31 @@ function select_numstate_best(E, sx, num_states) end end -function truncate_clustered_hamiltonian(cl_h, β, cs, result_folder, inst; tol=1e-6, iter=iter) +function truncate_clustered_hamiltonian( + cl_h, + β, + cs, + result_folder, + inst; + tol = 1e-6, + iter = iter, +) states = Dict() saved_states = load_file(joinpath(result_folder, "$(inst).jld2")) if saved_states == nothing new_cl_h = clustered_hamiltonian_2site(cl_h, β) - beliefs = belief_propagation(new_cl_h, β; tol=1e-6, iter=iter) - cl_h = truncate_clustered_hamiltonian_2site_BP(cl_h, beliefs, cs, result_folder, inst; beta=β) + beliefs = belief_propagation(new_cl_h, β; tol = 1e-6, iter = iter) + cl_h = truncate_clustered_hamiltonian_2site_BP( + cl_h, + beliefs, + cs, + result_folder, + inst; + beta = β, + ) else states = saved_states cl_h = truncate_clustered_hamiltonian(cl_h, states) end cl_h -end \ No newline at end of file +end diff --git a/src/utils.jl b/src/utils.jl index ba3f3f8..1a82549 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,7 +1,4 @@ -export - zephyr_to_linear, - unique_neighbors, - load_openGM +export zephyr_to_linear, unique_neighbors, load_openGM import Base.Prehashed using HDF5 @@ -14,7 +11,7 @@ Rewriten from Dwave-networkx m - Grid parameter for the Zephyr lattice. t - Tile parameter for the Zephyr lattice; must be even. """ -function zephyr_to_linear(m::Int, t::Int, q::NTuple{5, Int}) +function zephyr_to_linear(m::Int, t::Int, q::NTuple{5,Int}) M = 2 * m + 1 u, w, k, j, z = q (((u * M + w) * t + k) * 2 + j) * m + z + 1 @@ -104,7 +101,11 @@ Args: Returns: dictionary with factors and funcitons defining the energy functional. """ -function load_openGM(fname::String, Nx::Union{Integer, Nothing}=nothing, Ny::Union{Integer, Nothing}=nothing) +function load_openGM( + fname::String, + Nx::Union{Integer,Nothing} = nothing, + Ny::Union{Integer,Nothing} = nothing, +) file = h5open(fname, "r") file_keys = collect(keys(read(file))) @@ -129,13 +130,16 @@ function load_openGM(fname::String, Nx::Union{Integer, Nothing}=nothing, Ny::Uni nn = pop!(F) n = [] - for _ in 1:nn + for _ = 1:nn tt = pop!(F) ny, nx = divrem(tt, Nx) push!(n, ny, nx) end if length(n) == 4 - if abs(n[1] - n[3]) + abs(n[2] - n[4]) ∉ [1,2] || (abs(n[1] - n[3]) + abs(n[2] - n[4]) == 2 && (abs(n[1] - n[3]) == 2 || abs(n[2] - n[4] == 2))) + if abs(n[1] - n[3]) + abs(n[2] - n[4]) ∉ [1, 2] || ( + abs(n[1] - n[3]) + abs(n[2] - n[4]) == 2 && + (abs(n[1] - n[3]) == 2 || abs(n[2] - n[4] == 2)) + ) throw(ErrorException("Not nearest neighbour or diagonal neighbors")) end end @@ -163,28 +167,34 @@ function load_openGM(fname::String, Nx::Union{Integer, Nothing}=nothing, Ny::Uni nn = pop!(J) n = [] - for _ in 1:nn + for _ = 1:nn push!(n, pop!(J)) end upper = lower + prod(n) - functions[ii] = reshape(V[lower + 1:upper], reverse(n)...)' + functions[ii] = reshape(V[lower+1:upper], reverse(n)...)' lower = upper end - result = Dict("fun" => functions, "fac" => factors, "N" => reshape(N, (Ny, Nx)), "Nx" => Nx, "Ny" => Ny) + result = Dict( + "fun" => functions, + "fac" => factors, + "N" => reshape(N, (Ny, Nx)), + "Nx" => Nx, + "Ny" => Ny, + ) result end benchmark_names = Dict( - "penguin-small" => (240,320), - "palm-small" => (240,360), - "clownfish-small" => (240,360), - "crops-small" => (240,360), - "pfau-small" => (240,320), - "lake-small" => (240,360), - "snail" => (240,320), - "fourcolors" => (240,320), - "strawberry-glass-2-small" => (320,240) -) \ No newline at end of file + "penguin-small" => (240, 320), + "palm-small" => (240, 360), + "clownfish-small" => (240, 360), + "crops-small" => (240, 360), + "pfau-small" => (240, 320), + "lake-small" => (240, 360), + "snail" => (240, 320), + "fourcolors" => (240, 320), + "strawberry-glass-2-small" => (320, 240), +) diff --git a/test/bp_1site.jl b/test/bp_1site.jl index 5a679c5..c7f18d0 100644 --- a/test/bp_1site.jl +++ b/test/bp_1site.jl @@ -7,33 +7,28 @@ Instance below looks like this: """ function create_larger_example_clustered_hamiltonian_tree_basic() - instance = Dict( - (1, 1) => -0.50, - (2, 2) => 0.25, - (3, 3) => -0.30, - (4, 4) => 0.10, - (1, 2) => -0.23, - (1, 3) => 1.10, - (3, 4) => 0.71 - ) - - ig = ising_graph(instance) - - assignment_rule = Dict( - 1 => (1, 1, 1), - 2 => (1, 2, 1), - 3 => (2, 1, 1), - 4 => (2, 2, 2) - ) - - cl_h = clustered_hamiltonian( - ig, - Dict{NTuple{3, Int}, Int}(), - spectrum = full_spectrum, - cluster_assignment_rule = assignment_rule, - ) - - ig, cl_h + instance = Dict( + (1, 1) => -0.50, + (2, 2) => 0.25, + (3, 3) => -0.30, + (4, 4) => 0.10, + (1, 2) => -0.23, + (1, 3) => 1.10, + (3, 4) => 0.71, + ) + + ig = ising_graph(instance) + + assignment_rule = Dict(1 => (1, 1, 1), 2 => (1, 2, 1), 3 => (2, 1, 1), 4 => (2, 2, 2)) + + cl_h = clustered_hamiltonian( + ig, + Dict{NTuple{3,Int},Int}(), + spectrum = full_spectrum, + cluster_assignment_rule = assignment_rule, + ) + + ig, cl_h end """ @@ -46,47 +41,47 @@ Instance below looks like this: 7 -- 8 -- 9 """ function create_larger_example_clustered_hamiltonian_tree() - instance = Dict( - (1, 1) => 0.53, - (2, 2) => -0.25, - (3, 3) => 0.30, - (4, 4) => -0.10, - (5, 5) => -0.10, - (6, 6) => 0.10, - (8, 8) => 0.10, - (9, 9) => 0.01, - (1, 2) => -1.00, - (2, 3) => 1.00, - (1, 4) => 0.33, - (4, 5) => 0.76, - (5, 6) => -0.45, - (4, 7) => -0.28, - (7, 8) => 0.36, - (8, 9) => -1.07 - ) - - ig = ising_graph(instance) - - assignment_rule = Dict( - 1 => (1, 1, 1), - 2 => (1, 2, 1), - 3 => (1, 3, 1), - 4 => (2, 1, 1), - 5 => (2, 2, 1), - 6 => (2, 3, 1), - 7 => (3, 1, 1), - 8 => (3, 2, 1), - 9 => (3, 3, 1) - ) - - cl_h = clustered_hamiltonian( - ig, - Dict{NTuple{3, Int}, Int}(), - spectrum = full_spectrum, - cluster_assignment_rule = assignment_rule, - ) - - ig, cl_h + instance = Dict( + (1, 1) => 0.53, + (2, 2) => -0.25, + (3, 3) => 0.30, + (4, 4) => -0.10, + (5, 5) => -0.10, + (6, 6) => 0.10, + (8, 8) => 0.10, + (9, 9) => 0.01, + (1, 2) => -1.00, + (2, 3) => 1.00, + (1, 4) => 0.33, + (4, 5) => 0.76, + (5, 6) => -0.45, + (4, 7) => -0.28, + (7, 8) => 0.36, + (8, 9) => -1.07, + ) + + ig = ising_graph(instance) + + assignment_rule = Dict( + 1 => (1, 1, 1), + 2 => (1, 2, 1), + 3 => (1, 3, 1), + 4 => (2, 1, 1), + 5 => (2, 2, 1), + 6 => (2, 3, 1), + 7 => (3, 1, 1), + 8 => (3, 2, 1), + 9 => (3, 3, 1), + ) + + cl_h = clustered_hamiltonian( + ig, + Dict{NTuple{3,Int},Int}(), + spectrum = full_spectrum, + cluster_assignment_rule = assignment_rule, + ) + + ig, cl_h end """ @@ -98,63 +93,71 @@ Instance below looks like this: """ function create_larger_example_clustered_hamiltonian_tree_pathological() - instance = Dict( - (1, 1) => 0.52, - (2, 2) => 0.25, - (3, 3) => -0.31, - (4, 4) => 0.17, - (5, 5) => -0.12, - (6, 6) => 0.13, - (7, 7) => 0.00, - (8, 8) => 0.43, - (1, 2) => -1.01, - (1, 3) => 1.00, - (3, 4) => 0.97, - (3, 5) => -0.98, - (5, 6) => 1.00, - (2, 7) => 0.53, - (3, 7) => 1.06, - (5, 8) => -0.64, - ) - - ig = ising_graph(instance) - - assignment_rule = Dict( - 1 => (1, 1), - 2 => (1, 1), - 3 => (1, 1), - 4 => (1, 2), - 5 => (1, 2), - 6 => (1, 3), - 7 => (2, 1), - 8 => (2, 2) - ) - - cl_h = clustered_hamiltonian( - ig, - Dict{NTuple{2, Int}, Int}(), - spectrum = full_spectrum, - cluster_assignment_rule = assignment_rule, - ) - - ig, cl_h + instance = Dict( + (1, 1) => 0.52, + (2, 2) => 0.25, + (3, 3) => -0.31, + (4, 4) => 0.17, + (5, 5) => -0.12, + (6, 6) => 0.13, + (7, 7) => 0.00, + (8, 8) => 0.43, + (1, 2) => -1.01, + (1, 3) => 1.00, + (3, 4) => 0.97, + (3, 5) => -0.98, + (5, 6) => 1.00, + (2, 7) => 0.53, + (3, 7) => 1.06, + (5, 8) => -0.64, + ) + + ig = ising_graph(instance) + + assignment_rule = Dict( + 1 => (1, 1), + 2 => (1, 1), + 3 => (1, 1), + 4 => (1, 2), + 5 => (1, 2), + 6 => (1, 3), + 7 => (2, 1), + 8 => (2, 2), + ) + + cl_h = clustered_hamiltonian( + ig, + Dict{NTuple{2,Int},Int}(), + spectrum = full_spectrum, + cluster_assignment_rule = assignment_rule, + ) + + ig, cl_h end @testset "Belief propagation" begin - for (ig, cl_h) ∈ [create_larger_example_clustered_hamiltonian_tree_basic(), - create_larger_example_clustered_hamiltonian_tree(), - create_larger_example_clustered_hamiltonian_tree_pathological()] - for beta ∈ [0.5, 1] - iter = 16 - beliefs = belief_propagation(cl_h, beta; iter=iter) - exact_marginal = Dict() - for k in keys(beliefs) - push!(exact_marginal, k => [exact_cond_prob(cl_h, beta, Dict(k => a)) for a in 1:length(beliefs[k])]) - end - for v in keys(beliefs) - temp = -log.(exact_marginal[v])./beta - @test beliefs[v] ≈ temp .- minimum(temp) - end - end - end + for (ig, cl_h) ∈ [ + create_larger_example_clustered_hamiltonian_tree_basic(), + create_larger_example_clustered_hamiltonian_tree(), + create_larger_example_clustered_hamiltonian_tree_pathological(), + ] + for beta ∈ [0.5, 1] + iter = 16 + beliefs = belief_propagation(cl_h, beta; iter = iter) + exact_marginal = Dict() + for k in keys(beliefs) + push!( + exact_marginal, + k => [ + exact_cond_prob(cl_h, beta, Dict(k => a)) for + a = 1:length(beliefs[k]) + ], + ) + end + for v in keys(beliefs) + temp = -log.(exact_marginal[v]) ./ beta + @test beliefs[v] ≈ temp .- minimum(temp) + end + end + end end diff --git a/test/bp_2site.jl b/test/bp_2site.jl index 1d481f2..162f01b 100644 --- a/test/bp_2site.jl +++ b/test/bp_2site.jl @@ -9,63 +9,63 @@ Instance below looks like this: """ function create_larger_example_clustered_hamiltonian_tree_2site() - instance = Dict( - (1, 1) => 0.50, - (2, 2) => -0.25, - (3, 3) => 0.30, - (4, 4) => -0.10, - (5, 5) => 0.10, - (6, 6) => 0.20, - (7, 7) => -0.40, - (8, 8) => 0.50, - (1, 2) => -0.01, - (1, 3) => 1.00, - (3, 4) => 1.10, - (2, 5) => 0.77, - (6, 7) => -1.43, - (6, 8) => 0.21, - (6, 2) => 1.00 - ) - - ig = ising_graph(instance) - - assignment_rule1 = Dict( - 1 => (1, 1), - 2 => (1, 1), - 3 => (1, 2), - 4 => (1, 2), - 5 => (2, 1), - 6 => (2, 1), - 7 => (2, 2), - 8 => (2, 2) - ) - - assignment_rule2 = Dict( - 1 => (1, 1, 1), - 2 => (1, 1, 2), - 3 => (1, 2, 1), - 4 => (1, 2, 2), - 5 => (2, 1, 1), - 6 => (2, 1, 2), - 7 => (2, 2, 1), - 8 => (2, 2, 2) - ) - - cl_h1 = clustered_hamiltonian( - ig, - Dict{NTuple{2, Int}, Int}(), - spectrum = full_spectrum, - cluster_assignment_rule = assignment_rule1, - ) - - cl_h2 = clustered_hamiltonian( - ig, - Dict{NTuple{3, Int}, Int}(), - spectrum = full_spectrum, - cluster_assignment_rule = assignment_rule2, - ) - - ig, cl_h1, cl_h2 + instance = Dict( + (1, 1) => 0.50, + (2, 2) => -0.25, + (3, 3) => 0.30, + (4, 4) => -0.10, + (5, 5) => 0.10, + (6, 6) => 0.20, + (7, 7) => -0.40, + (8, 8) => 0.50, + (1, 2) => -0.01, + (1, 3) => 1.00, + (3, 4) => 1.10, + (2, 5) => 0.77, + (6, 7) => -1.43, + (6, 8) => 0.21, + (6, 2) => 1.00, + ) + + ig = ising_graph(instance) + + assignment_rule1 = Dict( + 1 => (1, 1), + 2 => (1, 1), + 3 => (1, 2), + 4 => (1, 2), + 5 => (2, 1), + 6 => (2, 1), + 7 => (2, 2), + 8 => (2, 2), + ) + + assignment_rule2 = Dict( + 1 => (1, 1, 1), + 2 => (1, 1, 2), + 3 => (1, 2, 1), + 4 => (1, 2, 2), + 5 => (2, 1, 1), + 6 => (2, 1, 2), + 7 => (2, 2, 1), + 8 => (2, 2, 2), + ) + + cl_h1 = clustered_hamiltonian( + ig, + Dict{NTuple{2,Int},Int}(), + spectrum = full_spectrum, + cluster_assignment_rule = assignment_rule1, + ) + + cl_h2 = clustered_hamiltonian( + ig, + Dict{NTuple{3,Int},Int}(), + spectrum = full_spectrum, + cluster_assignment_rule = assignment_rule2, + ) + + ig, cl_h1, cl_h2 end """ @@ -79,121 +79,130 @@ Instance below looks like this: """ function create_larger_example_clustered_hamiltonian_tree_2site_pathological() - instance = Dict( - (1, 1) => -0.50, - (2, 2) => 0.25, - (3, 3) => 0.31, - (4, 4) => 0.77, - (5, 5) => -0.15, - (6, 6) => 0.21, - (7, 7) => -0.10, - (8, 8) => 0.13, - (9, 9) => 0.01, - (10, 10) => 0.5, - (1, 2) => -0.01, - (1, 3) => 0.79, - (2, 5) => -0.93, - (5, 6) => 0.81, - (5, 9) => -0.17, - (3, 4) => 0.71, - (4, 6) => -0.43, - (5, 8) => 0.75, - (7, 10)=> -1.03, - (8, 9) => 0.65, - (9, 10)=> -0.57 - ) - - ig = ising_graph(instance) - - assignment_rule1 = Dict( - 1 => (1, 1), - 2 => (1, 1), - 3 => (1, 2), - 4 => (1, 2), - 5 => (1, 2), - 6 => (1, 3), - 7 => (2, 2), - 8 => (2, 2), - 9 => (2, 2), - 10 => (3, 2) - ) - - assignment_rule2 = Dict( - 1 => (1, 1, 1), - 2 => (1, 1, 2), - 3 => (1, 2, 1), - 4 => (1, 2, 1), - 5 => (1, 2, 2), - 6 => (1, 3, 2), - 7 => (2, 2, 1), - 8 => (2, 2, 1), - 9 => (2, 2, 2), - 10 => (3, 2, 2) - ) - - cl_h1 = clustered_hamiltonian( - ig, - Dict{NTuple{2, Int}, Int}(), - spectrum = full_spectrum, - cluster_assignment_rule = assignment_rule1, - ) - - cl_h2 = clustered_hamiltonian( - ig, - Dict{NTuple{3, Int}, Int}(), - spectrum = full_spectrum, - cluster_assignment_rule = assignment_rule2, - ) - - ig, cl_h1, cl_h2 + instance = Dict( + (1, 1) => -0.50, + (2, 2) => 0.25, + (3, 3) => 0.31, + (4, 4) => 0.77, + (5, 5) => -0.15, + (6, 6) => 0.21, + (7, 7) => -0.10, + (8, 8) => 0.13, + (9, 9) => 0.01, + (10, 10) => 0.5, + (1, 2) => -0.01, + (1, 3) => 0.79, + (2, 5) => -0.93, + (5, 6) => 0.81, + (5, 9) => -0.17, + (3, 4) => 0.71, + (4, 6) => -0.43, + (5, 8) => 0.75, + (7, 10) => -1.03, + (8, 9) => 0.65, + (9, 10) => -0.57, + ) + + ig = ising_graph(instance) + + assignment_rule1 = Dict( + 1 => (1, 1), + 2 => (1, 1), + 3 => (1, 2), + 4 => (1, 2), + 5 => (1, 2), + 6 => (1, 3), + 7 => (2, 2), + 8 => (2, 2), + 9 => (2, 2), + 10 => (3, 2), + ) + + assignment_rule2 = Dict( + 1 => (1, 1, 1), + 2 => (1, 1, 2), + 3 => (1, 2, 1), + 4 => (1, 2, 1), + 5 => (1, 2, 2), + 6 => (1, 3, 2), + 7 => (2, 2, 1), + 8 => (2, 2, 1), + 9 => (2, 2, 2), + 10 => (3, 2, 2), + ) + + cl_h1 = clustered_hamiltonian( + ig, + Dict{NTuple{2,Int},Int}(), + spectrum = full_spectrum, + cluster_assignment_rule = assignment_rule1, + ) + + cl_h2 = clustered_hamiltonian( + ig, + Dict{NTuple{3,Int},Int}(), + spectrum = full_spectrum, + cluster_assignment_rule = assignment_rule2, + ) + + ig, cl_h1, cl_h2 end @testset "Belief propagation 2site" begin - for (ig, cl_h1, cl_h2) ∈ [create_larger_example_clustered_hamiltonian_tree_2site(), - create_larger_example_clustered_hamiltonian_tree_2site_pathological()] - for beta ∈ [0.6, 1.1] - tol = 1e-12 - iter = 16 - num_states=10 - - new_cl_h1 = clustered_hamiltonian_2site(cl_h2, beta) - - @test vertices(new_cl_h1) == vertices(cl_h1) - @test edges(new_cl_h1) == edges(cl_h1) - for e ∈ vertices(new_cl_h1) - @test get_prop(new_cl_h1, e, :spectrum).energies ≈ get_prop(cl_h1, e, :spectrum).energies - end - for e ∈ edges(new_cl_h1) - E = get_prop(new_cl_h1, src(e), dst(e), :en) - @cast E[(l1, l2), (r1, r2)] := E.e11[l1, r1] + E.e21[l2, r1] + E.e12[l1, r2] + E.e22[l2, r2] - @test E == get_prop(cl_h1, src(e), dst(e), :en) - end - for e ∈ edges(new_cl_h1) - il1 = get_prop(new_cl_h1, src(e), dst(e), :ipl) - il2 = get_prop(cl_h1, src(e), dst(e), :ipl) - ir1 = get_prop(new_cl_h1, src(e), dst(e), :ipr) - ir2 = get_prop(cl_h1, src(e), dst(e), :ipr) - - pl1 = get_projector!(get_prop(new_cl_h1, :pool_of_projectors), il1, :CPU) - pl2 = get_projector!(get_prop(cl_h1, :pool_of_projectors), il2, :CPU) - pr1 = get_projector!(get_prop(new_cl_h1, :pool_of_projectors), ir1, :CPU) - pr2 = get_projector!(get_prop(cl_h1, :pool_of_projectors), ir2, :CPU) - @test pl1 == pl2 - @test pr1 == pr2 - end - - beliefs = belief_propagation(new_cl_h1, beta; iter=iter, tol=tol) - - exact_marginal = Dict() - for k in keys(beliefs) - temp = -1/beta .* log.([exact_cond_prob(cl_h1, beta, Dict(k => a)) for a in 1:length(beliefs[k])]) - push!(exact_marginal, k => temp .- minimum(temp)) - end - for v in keys(beliefs) - @test beliefs[v] ≈ exact_marginal[v] - end - end - end - end + for (ig, cl_h1, cl_h2) ∈ [ + create_larger_example_clustered_hamiltonian_tree_2site(), + create_larger_example_clustered_hamiltonian_tree_2site_pathological(), + ] + for beta ∈ [0.6, 1.1] + tol = 1e-12 + iter = 16 + num_states = 10 + + new_cl_h1 = clustered_hamiltonian_2site(cl_h2, beta) + + @test vertices(new_cl_h1) == vertices(cl_h1) + @test edges(new_cl_h1) == edges(cl_h1) + for e ∈ vertices(new_cl_h1) + @test get_prop(new_cl_h1, e, :spectrum).energies ≈ + get_prop(cl_h1, e, :spectrum).energies + end + for e ∈ edges(new_cl_h1) + E = get_prop(new_cl_h1, src(e), dst(e), :en) + @cast E[(l1, l2), (r1, r2)] := + E.e11[l1, r1] + E.e21[l2, r1] + E.e12[l1, r2] + E.e22[l2, r2] + @test E == get_prop(cl_h1, src(e), dst(e), :en) + end + for e ∈ edges(new_cl_h1) + il1 = get_prop(new_cl_h1, src(e), dst(e), :ipl) + il2 = get_prop(cl_h1, src(e), dst(e), :ipl) + ir1 = get_prop(new_cl_h1, src(e), dst(e), :ipr) + ir2 = get_prop(cl_h1, src(e), dst(e), :ipr) + + pl1 = get_projector!(get_prop(new_cl_h1, :pool_of_projectors), il1, :CPU) + pl2 = get_projector!(get_prop(cl_h1, :pool_of_projectors), il2, :CPU) + pr1 = get_projector!(get_prop(new_cl_h1, :pool_of_projectors), ir1, :CPU) + pr2 = get_projector!(get_prop(cl_h1, :pool_of_projectors), ir2, :CPU) + @test pl1 == pl2 + @test pr1 == pr2 + end + + beliefs = belief_propagation(new_cl_h1, beta; iter = iter, tol = tol) + + exact_marginal = Dict() + for k in keys(beliefs) + temp = + -1 / beta .* + log.([ + exact_cond_prob(cl_h1, beta, Dict(k => a)) for + a = 1:length(beliefs[k]) + ]) + push!(exact_marginal, k => temp .- minimum(temp)) + end + for v in keys(beliefs) + @test beliefs[v] ≈ exact_marginal[v] + end + end + end +end diff --git a/test/clustered_hamiltonian.jl b/test/clustered_hamiltonian.jl index 9389b6e..5bf9d25 100644 --- a/test/clustered_hamiltonian.jl +++ b/test/clustered_hamiltonian.jl @@ -5,28 +5,32 @@ using CSV enum(vec) = Dict(v => i for (i, v) ∈ enumerate(vec)) @testset "Lattice graph" begin - m = 4 - n = 4 - t = 4 - L = 128 + m = 4 + n = 4 + t = 4 + L = 128 - instance = "$(@__DIR__)/instances/chimera_droplets/$(L)power/001.txt" + instance = "$(@__DIR__)/instances/chimera_droplets/$(L)power/001.txt" - for T ∈ [Float16, Float32, Float64] + for T ∈ [Float16, Float32, Float64] ig = ising_graph(T, instance) - cl_h = clustered_hamiltonian(ig, 2, cluster_assignment_rule=super_square_lattice((m, n, 2*t))) + cl_h = clustered_hamiltonian( + ig, + 2, + cluster_assignment_rule = super_square_lattice((m, n, 2 * t)), + ) - @test collect(vertices(cl_h)) == [(i, j) for i ∈ 1:m for j ∈ 1:n] + @test collect(vertices(cl_h)) == [(i, j) for i ∈ 1:m for j ∈ 1:n] clv = [] cle = [] rank = rank_vec(ig) - for v ∈ vertices(cl_h) - cl = get_prop(cl_h, v, :cluster) - push!(clv, vertices(cl)) - push!(cle, collect(edges(cl))) + for v ∈ vertices(cl_h) + cl = get_prop(cl_h, v, :cluster) + push!(clv, vertices(cl)) + push!(cle, collect(edges(cl))) @test rank_vec(cl) == get_prop.(Ref(ig), vertices(cl), :rank) end @@ -38,58 +42,69 @@ enum(vec) = Dict(v => i for (i, v) ∈ enumerate(vec)) end @testset "Factor graph builds on pathological instance" begin - m = 3 - n = 4 - t = 3 - L = n * m * t - - instance = "$(@__DIR__)/instances/pathological/test_$(m)_$(n)_$(t).txt" - - ising = CSV.File(instance, types=[Int, Int, Float64], header=0, comment = "#") - - couplings = Dict((i, j) => v for (i, j, v) ∈ ising) - - cedges = Dict( - ((1, 1), (1, 2)) => [(1, 4), (1, 5), (1, 6)], - ((1, 1), (2, 1)) => [(1, 13)], - ((1, 2), (1, 3)) => [(4, 7), (5, 7), (6, 8), (6, 9)], - ((1, 2), (2, 2)) => [(6, 16), (6, 18), (5, 16)], - ((2, 1), (2, 2)) => [(13, 16), (13, 18)], - ((2, 2), (3, 2)) => [(18, 28)], - ((3, 2), (3, 3)) => [(28, 31), (28, 32), (28, 33), (29, 31), (29, 32), (29, 33), (30, 31), (30, 32), (30, 33)] - ) - - cells = Dict( - (1, 1) => [1], - (1, 2) => [4, 5, 6], - (1, 3) => [7, 8, 9], - (1, 4) => [], - (2, 1) => [13], - (2, 2) => [16, 18], - (2, 3) => [], - (2, 4) => [], - (3, 1) => [], - (3, 2) => [28, 29, 30], - (3, 3) => [31, 32, 33], - (3, 4) => [] - ) - - d = 2 - rank = Dict(c => fill(d, length(idx)) for (c, idx) ∈ cells if !isempty(idx)) - bond_dimensions = [2, 2, 8, 4, 2, 2, 8] - - for T ∈ [Float16, Float32, Float64] + m = 3 + n = 4 + t = 3 + L = n * m * t + + instance = "$(@__DIR__)/instances/pathological/test_$(m)_$(n)_$(t).txt" + + ising = CSV.File(instance, types = [Int, Int, Float64], header = 0, comment = "#") + + couplings = Dict((i, j) => v for (i, j, v) ∈ ising) + + cedges = Dict( + ((1, 1), (1, 2)) => [(1, 4), (1, 5), (1, 6)], + ((1, 1), (2, 1)) => [(1, 13)], + ((1, 2), (1, 3)) => [(4, 7), (5, 7), (6, 8), (6, 9)], + ((1, 2), (2, 2)) => [(6, 16), (6, 18), (5, 16)], + ((2, 1), (2, 2)) => [(13, 16), (13, 18)], + ((2, 2), (3, 2)) => [(18, 28)], + ((3, 2), (3, 3)) => [ + (28, 31), + (28, 32), + (28, 33), + (29, 31), + (29, 32), + (29, 33), + (30, 31), + (30, 32), + (30, 33), + ], + ) + + cells = Dict( + (1, 1) => [1], + (1, 2) => [4, 5, 6], + (1, 3) => [7, 8, 9], + (1, 4) => [], + (2, 1) => [13], + (2, 2) => [16, 18], + (2, 3) => [], + (2, 4) => [], + (3, 1) => [], + (3, 2) => [28, 29, 30], + (3, 3) => [31, 32, 33], + (3, 4) => [], + ) + + d = 2 + rank = Dict(c => fill(d, length(idx)) for (c, idx) ∈ cells if !isempty(idx)) + bond_dimensions = [2, 2, 8, 4, 2, 2, 8] + + for T ∈ [Float16, Float32, Float64] ig = ising_graph(T, instance) @test eltype(ig) == T cl_h = clustered_hamiltonian( ig, - spectrum=full_spectrum, - cluster_assignment_rule=super_square_lattice((m, n, t)), + spectrum = full_spectrum, + cluster_assignment_rule = super_square_lattice((m, n, t)), ) for (bd, e) in zip(bond_dimensions, edges(cl_h)) - ipl, en, ipr = get_prop(cl_h, e, :ipl), get_prop(cl_h, e, :en), get_prop(cl_h, e, :ipr) + ipl, en, ipr = + get_prop(cl_h, e, :ipl), get_prop(cl_h, e, :en), get_prop(cl_h, e, :ipr) pl = get_projector!(get_prop(cl_h, :pool_of_projectors), ipl, :CPU) pr = get_projector!(get_prop(cl_h, :pool_of_projectors), ipr, :CPU) @@ -99,7 +114,9 @@ end end for ((i, j), cedge) ∈ cedges - ipl, en, ipr = get_prop(cl_h, i, j, :ipl), get_prop(cl_h, i, j, :en), get_prop(cl_h, i, j, :ipr) + ipl, en, ipr = get_prop(cl_h, i, j, :ipl), + get_prop(cl_h, i, j, :en), + get_prop(cl_h, i, j, :ipr) pl = get_projector!(get_prop(cl_h, :pool_of_projectors), ipl, :CPU) pr = get_projector!(get_prop(cl_h, :pool_of_projectors), ipr, :CPU) base_i = all_states(rank[i]) @@ -135,15 +152,18 @@ end @testset "each edge comprises expected bunch of edges from source Ising graph" begin for e ∈ edges(cl_h) outer_edges = get_prop(cl_h, e, :outer_edges) - @test issetequal(cedges[(src(e), dst(e))], [(src(oe), dst(oe)) for oe ∈ outer_edges]) + @test issetequal( + cedges[(src(e), dst(e))], + [(src(oe), dst(oe)) for oe ∈ outer_edges], + ) end end end end -function create_example_clustered_hamiltonian(::Type{T}) where T +function create_example_clustered_hamiltonian(::Type{T}) where {T} J12 = -1 - h1 = 1/2 + h1 = 1 / 2 h2 = 0.75 D = Dict((1, 2) => J12, (1, 1) => h1, (2, 2) => h2) @@ -157,11 +177,12 @@ function create_example_clustered_hamiltonian(::Type{T}) where T ) end -cl_h_state_to_spin = [([1, 1], [-1, -1]), ([1, 2], [-1, 1]), ([2, 1], [1, -1]), ([2, 2], [1, 1])] +cl_h_state_to_spin = + [([1, 1], [-1, -1]), ([1, 2], [-1, 1]), ([2, 1], [1, -1]), ([2, 2], [1, 1])] @testset "Decoding solution gives correct spin assignment" begin - for T ∈ [Float16, Float32, Float64] + for T ∈ [Float16, Float32, Float64] cl_h = create_example_clustered_hamiltonian(T) @test all(eltype(get_prop(cl_h, e, :en)) == T for e ∈ edges(cl_h)) for (state, spin_values) ∈ cl_h_state_to_spin @@ -184,86 +205,90 @@ Instance below looks like this: And we group the following spins together: [1, 2, 4, 5], [3, 6], [7, 8], [9]. """ function create_larger_example_clustered_hamiltonian() - instance = Dict( - (1, 1) => 0.5, - (2, 2) => 0.25, - (3, 3) => 0.3, - (4, 4) => 0.1, - (5, 5) => 0.0, - (6, 6) => -2.0, - (7, 7) => -1.0, - (8, 8) => 2.0, - (9, 9) => 3.1, - (1, 2) => -1.0, - (2, 3) => 1.0, - (4, 5) => 0.5, - (5, 6) => -0.3, - (7, 8) => 0.1, - (8, 9) => 2.2, - (1, 4) => -1.7, - (4, 7) => 1.2, - (2, 5) => 0.2, - (5, 8) => 0.3, - (3, 6) => 1.1, - (6, 9) => 0.7 - ) - ig = ising_graph(instance) - - assignment_rule = Dict( - 1 => (1, 1), - 2 => (1, 1), - 4 => (1, 1), - 5 => (1, 1), - 3 => (1, 2), - 6 => (1, 2), - 7 => (2, 1), - 8 => (2, 1), - 9 => (2, 2) - ) - - cl_h = clustered_hamiltonian( - ig, - Dict{NTuple{2, Int}, Int}(), - spectrum = full_spectrum, - cluster_assignment_rule = assignment_rule, - ) - - ig, cl_h + instance = Dict( + (1, 1) => 0.5, + (2, 2) => 0.25, + (3, 3) => 0.3, + (4, 4) => 0.1, + (5, 5) => 0.0, + (6, 6) => -2.0, + (7, 7) => -1.0, + (8, 8) => 2.0, + (9, 9) => 3.1, + (1, 2) => -1.0, + (2, 3) => 1.0, + (4, 5) => 0.5, + (5, 6) => -0.3, + (7, 8) => 0.1, + (8, 9) => 2.2, + (1, 4) => -1.7, + (4, 7) => 1.2, + (2, 5) => 0.2, + (5, 8) => 0.3, + (3, 6) => 1.1, + (6, 9) => 0.7, + ) + ig = ising_graph(instance) + + assignment_rule = Dict( + 1 => (1, 1), + 2 => (1, 1), + 4 => (1, 1), + 5 => (1, 1), + 3 => (1, 2), + 6 => (1, 2), + 7 => (2, 1), + 8 => (2, 1), + 9 => (2, 2), + ) + + cl_h = clustered_hamiltonian( + ig, + Dict{NTuple{2,Int},Int}(), + spectrum = full_spectrum, + cluster_assignment_rule = assignment_rule, + ) + + ig, cl_h end function clustered_hamiltonian_energy(cl_h, state) - # This is highly inefficient, but simple, which makes it suitable for testing. - # If such a function is needed elsewhere, we need to implement it properly. - total_en = 0 - - # Collect local terms from each cluster - for (s, v) ∈ zip(state, vertices(cl_h)) - total_en += get_prop(cl_h, v, :spectrum).energies[s] - end - - # Collect inter-cluster terms - for edge ∈ edges(cl_h) - i, j = cl_h.reverse_label_map[src(edge)], cl_h.reverse_label_map[dst(edge)] - ipl, en, ipr = get_prop(cl_h, edge, :ipl), get_prop(cl_h, edge, :en), get_prop(cl_h, edge, :ipr) - pl = get_projector!(get_prop(cl_h, :pool_of_projectors), ipl, :CPU) - pr = get_projector!(get_prop(cl_h, :pool_of_projectors), ipr, :CPU) - edge_energy = en[pl, pr] - total_en += edge_energy[state[i], state[j]] - end - total_en + # This is highly inefficient, but simple, which makes it suitable for testing. + # If such a function is needed elsewhere, we need to implement it properly. + total_en = 0 + + # Collect local terms from each cluster + for (s, v) ∈ zip(state, vertices(cl_h)) + total_en += get_prop(cl_h, v, :spectrum).energies[s] + end + + # Collect inter-cluster terms + for edge ∈ edges(cl_h) + i, j = cl_h.reverse_label_map[src(edge)], cl_h.reverse_label_map[dst(edge)] + ipl, en, ipr = get_prop(cl_h, edge, :ipl), + get_prop(cl_h, edge, :en), + get_prop(cl_h, edge, :ipr) + pl = get_projector!(get_prop(cl_h, :pool_of_projectors), ipl, :CPU) + pr = get_projector!(get_prop(cl_h, :pool_of_projectors), ipr, :CPU) + edge_energy = en[pl, pr] + total_en += edge_energy[state[i], state[j]] + end + total_en end @testset "Decoding solution gives spins configuration with corresponding energies" begin - ig, cl_h = create_larger_example_clustered_hamiltonian() + ig, cl_h = create_larger_example_clustered_hamiltonian() - # Corresponding bases sizes for each cluster are 16, 4, 4, 2. - all_states = [[i, j, k, l] for i ∈ 1:16 for j ∈ 1:4 for k ∈ 1:4 for l ∈ 1:2] + # Corresponding bases sizes for each cluster are 16, 4, 4, 2. + all_states = [[i, j, k, l] for i ∈ 1:16 for j ∈ 1:4 for k ∈ 1:4 for l ∈ 1:2] - for state ∈ all_states - d = decode_clustered_hamiltonian_state(cl_h, state) - spins = zeros(length(d)) - for (k, v) ∈ d spins[k] = v end - σ = [Int.(spins)] - @test clustered_hamiltonian_energy(cl_h, state) ≈ energy(σ, ig)[] + for state ∈ all_states + d = decode_clustered_hamiltonian_state(cl_h, state) + spins = zeros(length(d)) + for (k, v) ∈ d + spins[k] = v + end + σ = [Int.(spins)] + @test clustered_hamiltonian_energy(cl_h, state) ≈ energy(σ, ig)[] end end diff --git a/test/ising.jl b/test/ising.jl index 439d273..763e280 100644 --- a/test/ising.jl +++ b/test/ising.jl @@ -6,7 +6,8 @@ using LabelledGraphs @testset "if input instance contains duplicate edges" begin for T ∈ [Float16, Float32, Float64] @test_throws ArgumentError ising_graph( - T, Dict((1, 1) => 2.0, (1, 2) => 0.5, (2, 1) => -1.0) + T, + Dict((1, 1) => 2.0, (1, 2) => 0.5, (2, 1) => -1.0), ) end end @@ -15,15 +16,27 @@ end for T ∈ [Float16, Float32, Float64] for (instance, source) ∈ ( ("$(@__DIR__)/instances/example.txt", "file"), - (Dict{Tuple{Int, Int}, T}((1, 1) => 0.1, (2, 2) => 0.5, (1, 4) => -2.0, (4, 2) => 1.0, (1, 2) => -0.3), "array") + ( + Dict{Tuple{Int,Int},T}( + (1, 1) => 0.1, + (2, 2) => 0.5, + (1, 4) => -2.0, + (4, 2) => 1.0, + (1, 2) => -0.3, + ), + "array", + ), ) @testset "Ising graph created from $(source)" begin expected_num_vertices = 3 - expected_biases = [T(1/10), T(1/2), T(0)] + expected_biases = [T(1 / 10), T(1 / 2), T(0)] expected_couplings = Dict( - LabelledEdge(1, 2) => -T(3/10), LabelledEdge(1, 4) => -T(2), LabelledEdge(2, 4) => T(1) + LabelledEdge(1, 2) => -T(3 / 10), + LabelledEdge(1, 4) => -T(2), + LabelledEdge(2, 4) => T(1), ) - expected_J_matrix = [[T(0) -T(3/10) -T(2)]; [T(0) T(0) T(1)]; [T(0) T(0) T(0)]] + expected_J_matrix = + [[T(0) -T(3 / 10) -T(2)]; [T(0) T(0) T(1)]; [T(0) T(0) T(0)]] ig = ising_graph(T, instance) @test eltype(ig) == T @@ -32,7 +45,8 @@ for T ∈ [Float16, Float32, Float64] end @testset "has collection of edges comprising all interactions from instance" begin # This test uses the fact that edges iterates in the lex ordering. - @test collect(edges(ig)) == [LabelledEdge(e...) for e ∈ [(1, 2), (1, 4), (2, 4)]] + @test collect(edges(ig)) == + [LabelledEdge(e...) for e ∈ [(1, 2), (1, 4), (2, 4)]] end @testset "stores biases as property of vertices" begin @test biases(ig) == expected_biases @@ -49,21 +63,17 @@ end @testset "Ising graph created with additional parameters" begin expected_biases = [-0.1, -0.5, 0.0] - expected_couplings = Dict( - Edge(1, 2) => 0.3, - Edge(1, 3) => 2.0, - Edge(2, 3) => -1.0 - ) + expected_couplings = Dict(Edge(1, 2) => 0.3, Edge(1, 3) => 2.0, Edge(2, 3) => -1.0) expected_J_matrix = [ - [0 0.3 2.0 ]; - [0 0 -1.0]; - [0 0 0]; + [0 0.3 2.0] + [0 0 -1.0] + [0 0 0] ] ig = ising_graph( "$(@__DIR__)/instances/example.txt", - scale=-1, - rank_override=Dict(1 => 3, 4 => 4) + scale = -1, + rank_override = Dict(1 => 3, 4 => 4), ) @testset "has rank overriden by rank_override dict" begin @@ -85,21 +95,27 @@ end ig = ising_graph(instance) @test nv(ig) == N - for i ∈ 1:N @test has_vertex(ig, i) end + for i ∈ 1:N + @test has_vertex(ig, i) + end A = adjacency_matrix(ig) B = zeros(Int, N, N) for i ∈ 1:N nbrs = SpinGlassNetworks.unique_neighbors(ig, i) - for j ∈ nbrs B[i, j] = 1 end + for j ∈ nbrs + B[i, j] = 1 + end end @test B + B' == A @testset "Reading from Dict" begin instance_dict = Dict() - ising = CSV.File(instance, types=[Int, Int, Float64], header=0, comment = "#") + ising = CSV.File(instance, types = [Int, Int, Float64], header = 0, comment = "#") - for (i, j, v) ∈ ising push!(instance_dict, (i, j) => v) end + for (i, j, v) ∈ ising + push!(instance_dict, (i, j) => v) + end ig = ising_graph(instance) ig_dict = ising_graph(instance_dict) @@ -116,18 +132,20 @@ end β = 1 instance = "$(@__DIR__)/instances/pathological/test_$(m)_$(n)_$(t).txt" - ising = CSV.File(instance, types=[Int, Int, Float64], header=0, comment = "#") + ising = CSV.File(instance, types = [Int, Int, Float64], header = 0, comment = "#") ig = ising_graph(instance) conf = [ [-1, 1, 1, -1, -1, -1, 1, 1, 1, -1, 1, 1, -1, 1, -1, 1], [-1, 1, 1, -1, -1, -1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1], [-1, 1, 1, -1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, 1], - [-1, 1, 1, -1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1] + [-1, 1, 1, -1, -1, 1, 1, 1, 1, -1, 1, 1, -1, 1, -1, -1], ] eng = energy(conf, ig) couplings = Dict() - for (i, j, v) ∈ ising push!(couplings, (i, j) => v) end + for (i, j, v) ∈ ising + push!(couplings, (i, j) => v) + end cedges = Dict() push!(cedges, (1, 2) => [(1, 4), (1, 5), (1, 6)]) @@ -138,9 +156,17 @@ end push!(cedges, (6, 10) => [(18, 28)]) push!( cedges, - (10, 11) => [(28, 31), (28, 32), (28, 33), (29, 31), (29, 32), - (29, 33),(30, 31), (30, 32), (30, 33) - ] + (10, 11) => [ + (28, 31), + (28, 32), + (28, 33), + (29, 31), + (29, 32), + (29, 33), + (30, 31), + (30, 32), + (30, 33), + ], ) push!(cedges, (2, 2) => [(4, 5), (4, 6), (5, 6), (6, 6)]) @@ -181,7 +207,7 @@ end push!(config, 30 => [-1, -1, -1, -1]) push!(config, 31 => [1, 1, 1, 1]) push!(config, 32 => [-1, -1, -1, -1]) - push!(config, 33 => [1,-1, 1, -1]) + push!(config, 33 => [1, -1, 1, -1]) push!(config, 34 => [0, 0, 0, 0]) push!(config, 35 => [0, 0, 0, 0]) push!(config, 36 => [0, 0, 0, 0]) @@ -189,13 +215,33 @@ end num_config = length(config[1]) exact_energy = _energy(config, couplings, cedges, num_config) - low_energies = [-16.4, -16.4, -16.4, -16.4, -16.1, - -16.1, -16.1, -16.1, -15.9, -15.9, - -15.9, -15.9, -15.9, -15.9, -15.6, - -15.6, -15.6, -15.6, -15.6, -15.6, - -15.4, -15.4 + low_energies = [ + -16.4, + -16.4, + -16.4, + -16.4, + -16.1, + -16.1, + -16.1, + -16.1, + -15.9, + -15.9, + -15.9, + -15.9, + -15.9, + -15.9, + -15.6, + -15.6, + -15.6, + -15.6, + -15.6, + -15.6, + -15.4, + -15.4, ] - for i ∈ 1:num_config @test exact_energy[i] == low_energies[i] == eng[i] end + for i ∈ 1:num_config + @test exact_energy[i] == low_energies[i] == eng[i] + end end end @@ -229,12 +275,12 @@ end @testset "Some vertices of degree zero, but nonzero field" begin instance = Dict( - (1, 1) => 0.1, - (2, 2) => 0.5, - (1, 4) => -2.0, - (4, 2) => 1.0, - (1, 2) => -0.3, - (5, 5) => 0.1 + (1, 1) => 0.1, + (2, 2) => 0.5, + (1, 4) => -2.0, + (4, 2) => 1.0, + (1, 2) => -0.3, + (5, 5) => 0.1, ) ig = ising_graph(instance) ng = prune(ig) @@ -244,12 +290,12 @@ end @testset "Some vertices of degree zero and zero field" begin instance = Dict( - (1, 1) => 0.1, - (2, 2) => 0.5, - (1, 4) => -2.0, - (4, 2) => 1.0, - (1, 2) => -0.3, - (5, 5) => 0.0 + (1, 1) => 0.1, + (2, 2) => 0.5, + (1, 4) => -2.0, + (4, 2) => 1.0, + (1, 2) => -0.3, + (5, 5) => 0.0, ) ig = ising_graph(instance) ng = prune(ig) diff --git a/test/runtests.jl b/test/runtests.jl index 4fe7616..37f7561 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -23,7 +23,7 @@ function _energy(config::Dict, couplings::Dict, cedges::Dict, n::Int) eng[m] += dot(s, J, r) end end - end + end end eng end @@ -34,8 +34,8 @@ my_tests = [ "bp_1site.jl", "bp_2site.jl", "utils.jl", - "spectrum.jl" - ] + "spectrum.jl", +] for my_test ∈ my_tests include(my_test) diff --git a/test/spectrum.jl b/test/spectrum.jl index 72992e9..cf7148a 100644 --- a/test/spectrum.jl +++ b/test/spectrum.jl @@ -10,7 +10,7 @@ using LabelledGraphs @testset "Naive brute force for +/-1" begin k = 2^N - sp = brute_force(ig, num_states=k) + sp = brute_force(ig, num_states = k) β = rand(Float64) ρ = gibbs_tensor(ig, β) @@ -37,7 +37,7 @@ using LabelledGraphs rank = get_prop(ig, :rank) all = prod(rank) - sp = full_spectrum(ig, num_states=all) + sp = full_spectrum(ig, num_states = all) β = rand(Float64) ρ = exp.(-β .* sp.energies) diff --git a/test/utils.jl b/test/utils.jl index b5b86b3..c822e8a 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -15,12 +15,12 @@ instances = ["P2"] #, "P4", "P8", "P16"] ig, max_cl_states, spectrum = brute_force, - cluster_assignment_rule = super_square_lattice((m, n, t)) + cluster_assignment_rule = super_square_lattice((m, n, t)), ) - @test nv(cl_h) == s ^ 2 + @test nv(cl_h) == s^2 if s > 1 - @test all(has_edge(cl_h, (l, k), (l+1, k-1)) for l ∈ 1:s-1, k ∈ 2:s) + @test all(has_edge(cl_h, (l, k), (l + 1, k - 1)) for l ∈ 1:s-1, k ∈ 2:s) end end end