Merge pull request #18 from lanl-ansi/fg-sampler

Generic Factor Graph Sampler

Author: marc-vuffray

src/GraphicalModelLearning.jl

@@ -146,7 +146,7 @@ function learn{T <: Real}(samples::Array{T,2}, formulation::multiRISE, method::N
         end
     end
 
-    return reconstruction
+    return FactorGraph(inter_order, num_spins, :spin, reconstruction) 
 end
 
 function learn{T <: Real}(samples::Array{T,2}, formulation::RISE, method::NLP)

src/models.jl

@@ -1,6 +1,6 @@
 # data structures graphical models
 
-export FactorGraph
+export FactorGraph, jsondata
 
 alphabets = [:spin, :boolean, :integer, :integer_pos, :real, :real_pos]
 
@@ -14,6 +14,9 @@ type FactorGraph{T <: Real}
 end
 FactorGraph{T <: Real}(order::Int, varible_count::Int, alphabet::Symbol, terms::Dict{Tuple,T}) = FactorGraph{T}(order, varible_count, alphabet, terms, Nullable{Vector{String}}())
 FactorGraph{T <: Real}(matrix::Array{T,2}) = convert(FactorGraph{T}, matrix)
+FactorGraph{T <: Real}(dict::Dict{Tuple,T}) = convert(FactorGraph{T}, dict)
+FactorGraph(list::Array{Any,1}) = convert(FactorGraph, list)
+
 
 function check_model_data{T <: Real}(order::Int, varible_count::Int, alphabet::Symbol, terms::Dict{Tuple,T}, variable_names::Nullable{Vector{String}})
     if !in(alphabet, alphabets)
@@ -35,11 +38,14 @@ function check_model_data{T <: Real}(order::Int, varible_count::Int, alphabet::S
                 error("a term has an index of $(index) but it should be in the range of 1:$(varible_count)")
                 return false
             end
+            #= 
+            # TODO see when this should be enforced
             if i > 1
                 if k[i-1] > index
                     error("the term $(k) does not have ascending indices")
                 end
             end
+            =#
         end
     end
     return true
@@ -53,12 +59,20 @@ function Base.show(io::IO, gm::FactorGraph)
         println(io, "  ", get(gm.variable_names))
     end
 
-    println(io, "terms: ")
-    for k in sort(collect(keys(gm.terms)))
+    println(io, "terms: $(length(gm.terms))")
+    for k in sort(collect(keys(gm.terms)), by=(x)->(length(x),x))
         println("  ", k, " => ", gm.terms[k])
     end
 end
 
+function jsondata{T <: Real}(gm::FactorGraph{T})
+    data = []
+    for k in sort(collect(keys(gm.terms)), by=(x)->(length(x),x))
+        push!(data, Dict("term" => k, "weight" => gm.terms[k]))
+    end
+    return data
+end
+
 Base.start(gm::FactorGraph) = start(gm.terms)
 Base.next(gm::FactorGraph, state) = next(gm.terms, state)
 Base.done(gm::FactorGraph, state) = done(gm.terms, state)
@@ -68,7 +82,6 @@ Base.length(gm::FactorGraph) = length(gm.terms)
 Base.getindex(gm::FactorGraph, i) = gm.terms[i]
 Base.keys(gm::FactorGraph) = keys(gm.terms)
 
-
 function diag_keys(gm::FactorGraph)
     dkeys = Tuple[]
     for i in 1:gm.varible_count
@@ -167,6 +180,48 @@ function Base.convert{T <: Real}(::Type{Dict{Tuple,T}}, m::Array{T,2})
 end
 
 
+function Base.convert(::Type{FactorGraph}, list::Array{Any,1})
+    info("assuming spin alphabet")
+    alphabet = :spin
+
+    max_variable = 0
+    max_order = 0
+    terms = Dict{Tuple,Float64}()
+
+    for item in list
+        term = item["term"]
+        weight = item["weight"]
+        terms[tuple(term...)] = weight
+
+        @assert minimum(term) > 0
+        max_order = max(max_order, length(term))
+        max_variable = max(max_variable, maximum(term))
+    end
+
+    info("dectected $(max_variable) variables with order $(max_order)")
+
+    return FactorGraph(max_order, max_variable, alphabet, terms)
+end
+
+
+Base.convert{T <: Real}(::Type{FactorGraph}, dict::Dict{Tuple,T}) = convert(FactorGraph{T}, dict)
+function Base.convert{T <: Real}(::Type{FactorGraph{T}}, dict::Dict{Tuple,T})
+    info("assuming spin alphabet")
+    alphabet = :spin
+
+    max_variable = 0
+    max_order = 0
+    for (term,weight) in dict
+        @assert minimum(term) > 0
+        max_order = max(max_order, length(term))
+        max_variable = max(max_variable, maximum(term))
+    end
+
+    info("dectected $(max_variable) variables with order $(max_order)")
+
+    return FactorGraph(max_order, max_variable, alphabet, dict)
+end
+
 
 permutations(items, order::Int; asymmetric::Bool = false) = sort(permutations([], items, order, asymmetric))
 

src/sampling.jl

@@ -30,7 +30,8 @@ function weigh_proba{T <: Real}(int_representation::Int, adj::Array{T,2}, prior:
 end
 
 
-function sample_generation{T <: Real}(gm::FactorGraph{T}, samples_per_bin::Integer, bins::Int)
+# assumes second order
+function sample_generation_ising{T <: Real}(gm::FactorGraph{T}, samples_per_bin::Integer, bins::Int)
     @assert bins >= 1
 
     spin_number   = gm.varible_count
@@ -43,7 +44,37 @@ function sample_generation{T <: Real}(gm::FactorGraph{T}, samples_per_bin::Integ
     assignment_tmp = [0 for i in 1:spin_number] # pre allocate assignment memory
     weights = [weigh_proba(i, adjacency_matrix, prior_vector, assignment_tmp) for i in (0:config_number-1)]
 
-    items = [i for i in 0:(config_number-1)]
+    raw_sample = StatsBase.sample(items, StatsBase.Weights(weights), samples_per_bin*bins, ordered=false)
+    raw_sample_bins = reshape(raw_sample, bins, samples_per_bin)
+
+    spin_samples = []
+    for b in 1:bins
+        raw_binning = countmap(raw_sample_bins[b,:])
+        spin_sample = [ vcat(raw_binning[i], int_to_spin(i, spin_number)) for i in keys(raw_binning)]
+        push!(spin_samples, hcat(spin_sample...)')
+    end
+    return spin_samples
+end
+
+
+function weigh_proba{T <: Real}(int_representation::Int, gm::FactorGraph{T}, spins::Array{Int,1})
+    digits!(spins, int_representation, 2)
+    spins .= bool_to_spin.(spins)
+    evaluation = sum( weight*prod(spins[i] for i in term) for (term, weight) in gm) 
+    return exp(evaluation)
+end
+
+function sample_generation{T <: Real}(gm::FactorGraph{T}, samples_per_bin::Integer, bins::Int)
+    @assert bins >= 1
+    #info("use general sample model")
+
+    spin_number   = gm.varible_count
+    config_number = 2^spin_number
+
+    items   = [i for i in 0:(config_number-1)]
+    assignment_tmp = [0 for i in 1:spin_number] # pre allocate assignment memory
+    weights = [weigh_proba(i, gm, assignment_tmp) for i in (0:config_number-1)]
+
     raw_sample = StatsBase.sample(items, StatsBase.Weights(weights), samples_per_bin*bins, ordered=false)
     raw_sample_bins = reshape(raw_sample, bins, samples_per_bin)
 
@@ -61,14 +92,15 @@ sample{T <: Real}(gm::FactorGraph{T}, number_sample::Integer, replicates::Intege
 
 
 function sample{T <: Real}(gm::FactorGraph{T}, number_sample::Integer, replicates::Integer, sampler::Gibbs)
-    if gm.order != 2
-        error("sampling is only supported for FactorGraphs of order 2, given order $(gm.order)")
-    end
     if gm.alphabet != :spin
         error("sampling is only supported for spin FactorGraphs, given alphabet $(gm.alphabet)")
     end
 
-    samples = sample_generation(gm, number_sample, replicates)
+    if gm.order <= 2
+        samples = sample_generation_ising(gm, number_sample, replicates)
+    else
+        samples = sample_generation(gm, number_sample, replicates)
+    end
 
     return samples
 end

test/runtests.jl

@@ -22,6 +22,17 @@ end
 
 
 @testset "gibbs sampler" begin
+    for (name, gm) in gms
+        gm_tmp = deepcopy(gm)
+        gm_tmp.order = 3
+        srand(0) # fix random number generator
+        samples = sample(gm_tmp, gibbs_test_samples)
+        base_samples = readcsv("data/$(name)_samples.csv")
+        @test isapprox(samples, base_samples)
+    end
+end
+
+@testset "gibbs sampler, 2nd order" begin
     for (name, gm) in gms
         srand(0) # fix random number generator
         samples = sample(gm, gibbs_test_samples)