Package ising code

Author: guilgautier

src/PartialRejectionSampling.jl

@@ -1,5 +1,19 @@
 module PartialRejectionSampling
 
-# Write your package code here.
+using LinearAlgebra
+const LA = LinearAlgebra
+
+using Random, Distributions
+
+using LightGraphs
+const LG = LightGraphs
+
+using SimpleWeightedGraphs
+const SWG = SimpleWeightedGraphs
+
+using Plots, GraphPlot, Colors
+
+include("ising.jl")
+include("utils.jl")
 
 end

src/ising.jl

@@ -1,93 +1,182 @@
-import Random
-getRNG(seed::Integer = -1) = seed >= 0 ? Random.MersenneTwister(seed) : Random.GLOBAL_RNG
-getRNG(seed::Union{Random.MersenneTwister,Random._GLOBAL_RNG}) = seed
+struct Ising
+    "Grid dimension"
+    dims::AbstractVector{Int}
+    "Interaction graph"
+    g::LG.SimpleGraph{Int}
+    "Magnetization"
+    h::Union{Real, AbstractVector{Real}}
+    "Interaction"
+    J::Real
+end
 
+function Ising(
+    dims::AbstractVector{Int},
+    periodic::Bool,
+    h::Union{Real, AbstractVector{Real}},
+    J::Real
+)
+    if isa(h, AbstractVector)
+        @assert length(h) != prod(dims) "length(h) != prod(dims)"
+    end
+    g = LG.grid(dims; periodic=periodic)
+    return Ising(dims, g, h, J)
+end
 
-import LightGraphs
-const LG = LightGraphs
+function Ising(
+    dims::AbstractVector{Int},
+    periodic::Bool,
+    J::Real
+)
+    return Ising(dims, periodic, 0.0, J)
+end
 
-import SimpleWeightedGraphs
-const SWG = SimpleWeightedGraphs
+function energy(
+    ising::Ising,
+    state::AbstractVector{Int}
+)::Real
+    @assert all(x -> x ∈ [-1, 1], state) "state values not all {-1, 1}"
 
-using GraphPlot, Colors
+    E = 0.0
 
-function ising_dependency_graph(
-    dims::Vector{T};
-    periodic::Bool = false,
-    rng = -1,
-)::SWG.SimpleWeightedGraph{Int64,Float64} where {T}
-    rng = getRNG(rng)
-    g = SWG.SimpleWeightedGraph(LG.grid(dims, periodic = periodic))
-    for e in LG.edges(g)
+    h = ising.h
+    if isa(h, Real)
+        E = - h * sum(state)
+    else  # isa(h, AbstractVector)
+        E = - h' * state
+    end
+
+    J = ising.J
+    for e in LG.edges(ising.g)
         i, j = Tuple(e)
-        @inbounds g.weights[i, j] = g.weights[j, i] = rand(rng)
+        xᵢ, xⱼ = state[i], state[j]
+        E += J * xᵢ * xⱼ
     end
-    return g
+    return E
 end
 
-"""
-Find bad
-U_ij > exp(- |J| (1 - sign(J) x_i x_j))
-log(U_ij) > - |J| (1 - sign(J) x_i x_j)
-sign(J) x_i x_j < 0 & log(U_ij) > - 2 |J|
-"""
-function ising_find_bad_states(
-    g::SWG.SimpleWeightedGraph{T, U},
-    states::Vector{T},
-    J::U;
-    rng = -1,
-)::Set{T} where {T, U}
+function plot(
+    ising::Ising,
+    state::AbstractVector{Int}
+)
+
+    pos = collect(Iterators.product(1:ising.dims[1], 1:ising.dims[2]))[:]
+    locs_x, locs_y = map(x->x[1], pos), map(x->x[2], pos)
+
+    col_nodes = ifelse.(
+                    state .== 1,
+                    Colors.colorant"gray",
+                    Colors.colorant"white")
+
+    p = GraphPlot.gplot(
+            ising.g,
+            locs_x,
+            reverse(locs_y),
+            nodefillc=col_nodes)
+    #     nodelabel=LG.vertices(g),
+    #     arrowlengthfrac=0.05
+    #     edgestrokec=col_edges
+
+    display(p)
+end
+
+function sample_state!(
+    state::AbstractVector{T},
+    ising::Ising,
+    i::T;
+    rng=-1
+) where {T<:Int}
+    rng = getRNG(rng)
+    hᵢ = isa(ising.h, Real) ? ising.h : ising.h[i]
+    state[i] = rand(rng) < sigmoid(hᵢ) ? 1 : -1
+end
+
+function sample_states!(
+    state::AbstractVector{T},
+    ising::Ising,
+    indices::Set{T};
+    rng=-1
+) where {T<:Int}
+    for i in indices
+        sample_state!(state, ising, i; rng=rng)
+    end
+end
+
+function sample_conditional!(
+    state::AbstractVector{T},
+    ising::Ising,
+    i::T;
+    rng=-1
+) where {T<:Int}
+    rng = getRNG(rng)
+    hᵢ = isa(ising.h, Real) ? ising.h : ising.h[i]
+    proba = sigmoid(hᵢ + ising.J * sum(state[j] for j in LG.neighbors(ising.g, i)))
+    state[i] = rand(rng) < proba ? 1 : -1
+end
 
-    @assert LG.nv(g) == length(states)
+## Partial Rejection Sampling sampler (mix of Guo Jerrum Liu and Moka Kroese)
 
+function weighted_graph(
+    ising::Ising;
+    rng=-1
+)::SWG.SimpleWeightedGraph{Int64,Float64}
     rng = getRNG(rng)
+    wg = SWG.SimpleWeightedGraph(ising.g)
+    for e in LG.edges(wg)
+        i, j = Tuple(e)
+        @inbounds wg.weights[i, j] = wg.weights[j, i] = rand(rng)
+    end
+    return wg
+end
+
+function bad_states(
+    ising::Ising,
+    state::AbstractVector{T},
+    wg::SWG.SimpleWeightedGraph{T, U}
+)::Set{T} where {T, U}
 
-    _2absJ, sign_J = -2.0 * abs(J), sign(J)
+    _2absJ, sign_J = -2.0 * abs(ising.J), sign(ising.J)
 
     bad = Set{T}()
-    for e ∈ LG.edges(g)
+    for e ∈ LG.edges(wg)
         i, j = Tuple(e)
         # Break constraint:
         # log(U_ij) > - |J| (1 - sign(J) x_i x_j)
         # <=> sign(J) x_i x_j < 0 & log(U_ij) > - 2 |J|
-        if ((sign_J * states[i] * states[j]) < 0.0) & (log(e.weight) > _2absJ)
+        if ((sign_J * state[i] * state[j]) < 0.0) & (log(e.weight) > _2absJ)
             union!(bad, i, j)
         end
     end
     return bad
 end
 
-"""
-Find Res
-"""
-function ising_find_states_to_resample!(
-    g::SWG.SimpleWeightedGraph{T, U},
-    states::Vector{T},
-    J::U;
-    rng = -1
+function resampling_states!(
+    wg::SWG.SimpleWeightedGraph{T, U},
+    ising::Ising,
+    state::AbstractVector{T};
+    rng=-1
 )::Set{T} where {T, U}
 
     rng = getRNG(rng)
-    sign_J = sign(J)
-    R = ising_find_bad_states(g, states, J, rng = rng)
+    sign_J = sign(ising.J)
+    R = bad_states(ising, state, wg)
     ∂R, ∂R_tmp = copy(R), Set{T}()
     while !isempty(∂R)
         for i ∈ ∂R
-            for j ∈ LG.neighbors(g, i)
+            for j ∈ LG.neighbors(wg, i)
                 # Break constraint:
                 # log(U_ij) > - |J| (1 - sign(J) x_i x_j)
                 # <=> sign(J) x_i x_j < 0 & log(U_ij) > - 2 |J|
                 if j ∈ R
-                    if (sign_J * states[i] * states[j]) < 0.0
+                    if (sign_J * state[i] * state[j]) < 0.0
                         # U_ij can be increased
-                        @inbounds g.weights[i, j] = g.weights[j, i] = rand(rng)
+                        @inbounds wg.weights[i, j] = wg.weights[j, i] = rand(rng)
                     end
                 else # if j ∉ R
                     # x_j can be flipped to make sign(J) x_i x_j < 0
                     union!(R, j)
                     union!(∂R_tmp, j)
                     # followed by an increase of U_ij
-                    @inbounds g.weights[i, j] = g.weights[j, i] = rand(rng)
+                    @inbounds wg.weights[i, j] = wg.weights[j, i] = rand(rng)
                 end
             end
         end
@@ -96,85 +185,67 @@ function ising_find_states_to_resample!(
     return R
 end
 
-"""
-Resample
-"""
-function ising_sample_states!(
-    states::Vector{T},
-    res_ind::Set{T},
-    probas::Vector{Float64};
-    rng = -1,
-) where {T}
-    n = length(states)
-    @assert n == length(probas)
-    rng = getRNG(rng)
-    for i in res_ind
-        states[i] = rand(rng) < probas[i] ? 1 : -1
-    end
-end
-
-sigmoid(x) = @. 1 / (1 + exp(-x))
-
-function ising_prs(
-    dims::Vector{T},
-    h::Vector{U},
-    J::U;
-    periodic::Bool=false,
+function prs(
+    ising::Ising;
     rng=-1,
-) where {T, U}
+)::Tuple{AbstractVector{Int}, Int}
     rng = getRNG(rng)
 
-    g = ising_dependency_graph(dims, periodic=periodic, rng=rng)
-
-    n = LG.nv(g)
-    states = Vector{Int64}(undef, n)
-    probas = sigmoid.(2.0 .* h)
-    res = Set{T}(1:LG.nv(g))
+    g = weighted_graph(ising; rng=rng)
+    state = Vector{Int}(undef, LG.nv(g))
+    res = Set{Int}(1:LG.nv(g))
 
     cnt = 0
     while !isempty(res)
-        ising_sample_states!(states, res, probas, rng=rng)
-        res = ising_find_states_to_resample!(g, states, J, rng = rng)
+        sample_states!(state, ising, res; rng=rng)
+        res = resampling_states!(g, ising, state; rng=rng)
         cnt += 1
     end
 
-    return g, states, cnt
+    return state, cnt
 end
 
-function ising_prs(
-    dims::Vector{T},
-    h::U,
-    J::U;
-    periodic::Bool = false,
-    rng = -1,
-) where {T, U}
-    h_vec = fill(h, prod(dims))
-    return ising_prs(dims, h_vec, J, periodic=periodic, rng=rng)
-end
+## Perfect Gibbs sampler (Feng, Guo, Yin) https://arxiv.org/pdf/1907.06033.pdf
 
-function plot_ising(g, dims, state)
-    pos = collect(Iterators.product(1:dims[1], 1:dims[2]))[:]
-    locs_x, locs_y = map(x->x[1], pos), map(x->x[2], pos)
+function bayes_filter(
+    ising::Ising,
+    state::AbstractVector{T},
+    i::T,
+    R::Set{T};
+    rng=-1
+)::Bool where {T<:Int}
 
-    col_nodes = ifelse.(state .== 1, colorant"gray", colorant"white")
+    ∂i_R̄ = setdiff(LG.neighbors(ising.g, i), R)
+    isempty(∂i_R̄) && return true
 
-    p = gplot(g,
-        locs_x, reverse(locs_y),
-        nodefillc=col_nodes,
-    #     nodelabel=LG.vertices(g),
-    #     arrowlengthfrac=0.05
-    #     edgestrokec=col_edges
-    )
-    display(p)
+    rng = getRNG(rng)
+    Xi_J = state[i] * ising.J
+    acc_r = Xi_J * (-sign(Xi_J) * length(∂i_R̄) - sum(state[j] for j in ∂i_R̄))
+    return log(rand(rng)) < acc_r
 end
 
-dims = [14, 14] # if > (14, 14) the display becomes all black, don't know why !
-H, J = 0.0, -0.02 # Use Float
+function gibbs_perfect(
+    ising::Ising;
+    rng=-1
+)::Tuple{AbstractVector{Int}, Int}
+    rng = getRNG(rng)
 
-periodic = false
-seed = -1
-g, config, cnt = ising_prs(dims, H, J; periodic=periodic, rng=seed)
+    n = LG.nv(ising.g)
+    R = Set{Int}(1:n)
+    state = Vector{Int}(undef, n)
+    sample_states!(state, ising, R; rng=rng)
+
+    cnt = 1
+    while !isempty(R)
+        i = rand(rng, R)
+    if bayes_filter(ising, state, i, R; rng=-1)
+            sample_conditional!(state, ising, i; rng=rng)
+            delete!(R, i)
+        else
+            union!(R, LG.neighbors(ising.g, i))
+        end
+        cnt += 1
+    end
 
-plot_ising(g, dims, config)
-println("Number of resampling steps")
-print(cnt)
+    return state, cnt
+end

src/utils.jl

@@ -0,0 +1,4 @@
+getRNG(seed::Integer = -1) = seed >= 0 ? Random.MersenneTwister(seed) : Random.GLOBAL_RNG
+getRNG(seed::Union{Random.MersenneTwister,Random._GLOBAL_RNG}) = seed
+
+sigmoid(x) = @. 1 / (1 + exp(-x))