Adding Support for Factor Graphs (#13)

* adding basic factor graph data structure
* adding tests for factor graphs

Author: ccoffrin

README.md

@@ -13,11 +13,11 @@ Try the following commands in julia,
 ```
 using GraphicalModelLearning
 
-model = [0.0 0.1 0.2; 0.1 0.0 0.3; 0.2 0.3 0.0]
+model = FactorGraph([0.0 0.1 0.2; 0.1 0.0 0.3; 0.2 0.3 0.0])
 samples = sample(model, 100000)
 learned = learn(samples)
 
-err = abs.(model - learned)
+err = abs.(convert(Array{Float64,2}, model) - learned)
 ```
 
 Note that the first invocation of `learn` will be slow as the dependent libraries are compiled.  Subsequent calls will be fast.

src/GraphicalModelLearning.jl

@@ -13,6 +13,8 @@ using Ipopt
 
 using Compat # used for julia v0.5 abstract types
 
+include("models.jl")
+
 include("sampling.jl")
 
 @compat abstract type GMLFormulation end

src/models.jl

@@ -0,0 +1,150 @@
+# data structures graphical models
+
+export FactorGraph
+
+alphabets = [:spin, :boolean, :integer, :integer_pos, :real, :real_pos]
+
+type FactorGraph{T <: Real}
+    order::Int
+    varible_count::Int
+    alphabet::Symbol
+    terms::Dict{Tuple,T} # TODO, would be nice to have a stronger tuple type here
+    variable_names::Nullable{Vector{String}}
+    FactorGraph(a,b,c,d,e) = check_model_data(a,b,c,d,e) ? new(a,b,c,d,e) : error("generic init problem")
+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)
+
+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)
+        error("alphabet $(alphabet) is not supported")
+        return false 
+    end
+    if !isnull(variable_names) && length(variable_names) != varible_count
+        error("expected $(varible_count) but only given $(length(variable_names))")
+        return false 
+    end
+    for (k,v) in terms
+        if length(k) != order
+            error("a term has $(length(k)) indices but should have $(order) indices")
+            return false
+        end
+        for (i,index) in enumerate(k)
+            #println(i," ",index)
+            if index < 1 || index > varible_count
+                error("a term has an index of $(index) but it should be in the range of 1:$(varible_count)")
+                return false
+            end
+            if i > 1
+                if k[i-1] > index
+                    error("the term $(k) does not have ascending indices")
+                end
+            end
+        end
+    end
+    return true
+end
+
+function Base.show(io::IO, gm::FactorGraph)
+    println(io, "alphabet: ", gm.alphabet)
+    println(io, "vars: ", gm.varible_count)
+    if !isnull(gm.variable_names)
+        println(io, "variable names: ")
+        println(io, "  ", get(gm.variable_names))
+    end
+
+    println(io, "terms: ")
+    for k in sort(collect(keys(gm.terms)))
+        println("  ", k, " => ", gm.terms[k])
+    end
+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)
+
+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
+        key = diag_key(gm, i)
+        if key in keys(gm.terms)
+            push!(dkeys, key)
+        end
+    end
+    return sort(dkeys)
+end
+
+diag_key(gm::FactorGraph, i::Int) = tuple(fill(i, gm.order)...)
+
+#Base.diag{T <: Real}(gm::FactorGraph{T}) = [ get(gm.terms, diag_key(gm, i), zero(T)) for i in 1:gm.varible_count ]
+
+Base.DataFmt.writecsv{T <: Real}(io, gm::FactorGraph{T}, args...; kwargs...) = writecsv(io, convert(Array{T,2}, gm), args...; kwargs...)
+
+Base.convert{T <: Real}(::Type{FactorGraph}, m::Array{T,2}) = convert(FactorGraph{T}, m)
+function Base.convert{T <: Real}(::Type{FactorGraph{T}}, m::Array{T,2})
+    @assert size(m,1) == size(m,2) #check matrix is square
+
+    info("assuming spin alphabet")
+    alphabet = :spin
+    varible_count = size(m,1)
+
+    terms = Dict{Tuple,T}()
+    for key in permutations(1:varible_count, 2)
+        weight = m[key...]
+        if !isapprox(weight, 0.0)
+            terms[key] = weight
+        end
+
+        rev = reverse(key)
+        if !isapprox(m[rev...], 0.0) && !isapprox(m[key...], m[rev...])
+            delta = abs(m[key...] - m[rev...])
+            warn("values at $(key) and $(rev) differ by $(delta), only $(key) will be used")
+        end
+    end
+
+    return FactorGraph(2, varible_count, alphabet, terms)
+end
+
+function Base.convert{T <: Real}(::Type{Array{T,2}}, gm::FactorGraph{T})
+    if gm.order != 2
+        error("cannot convert a FactorGraph of order $(gm.order) to a matrix")
+    end
+
+    matrix = zeros(gm.varible_count, gm.varible_count)
+    for (k,v) in gm
+        matrix[k...] = v
+        r = reverse(k)
+        matrix[r...] = v
+    end
+
+    return matrix
+end
+
+
+
+permutations(items, order::Int; asymmetric::Bool = false) = sort(permutations([], items, order, asymmetric))
+
+function permutations(partical_perm::Array{Any,1}, items, order::Int, asymmetric::Bool)
+    if order == 0
+        return [tuple(partical_perm...)]
+    else
+        perms = []
+        for item in items
+            if !asymmetric && length(partical_perm) > 0 
+                if partical_perm[end] < item
+                    continue
+                end
+            end
+            perm = permutations(vcat([item], partical_perm), items, order-1, asymmetric)
+            append!(perms, perm)
+        end
+        return perms
+    end
+end
+

src/sampling.jl

@@ -13,30 +13,35 @@ function int_to_spin(int_representation::Int, spin_number::Int)
     return spin
 end
 
+
 function weigh_proba{T <: Real}(int_representation::Int, adj::Array{T,2}, prior::Array{T,1})
-    spin_number  = size(adj,1)
+    spin_number = size(adj,1)
     spins = int_to_spin(int_representation, spin_number)
     return exp(((0.5) * spins' * adj * spins + prior' * spins)[1])
 end
 
 
 bool_to_spin(bool::Int) = 2*bool-1
 
-function weigh_proba{T <: Real}(int_representation::Int, adj::Array{T,2}, prior::Array{T,1}, assignment_tmp::Array{Int,1})
-    digits!(assignment_tmp, int_representation, 2)
-    assignment_tmp .= bool_to_spin.(assignment_tmp)
-    return exp(((0.5) * assignment_tmp' * adj * assignment_tmp + prior' * assignment_tmp)[1])
+function weigh_proba{T <: Real}(int_representation::Int, adj::Array{T,2}, prior::Array{T,1}, spins::Array{Int,1})
+    digits!(spins, int_representation, 2)
+    spins .= bool_to_spin.(spins)
+    return exp(((0.5) * spins' * adj * spins + prior' * spins)[1])
 end
 
 
-function sample_generation{T <: Real}(samples_per_bin::Integer, adj::Array{T,2}, prior::Array{T,1}, bins::Int)
+function sample_generation{T <: Real}(gm::FactorGraph{T}, samples_per_bin::Integer, bins::Int)
     @assert bins >= 1
 
-    spin_number   = size(adj,1)
+    spin_number   = gm.varible_count
     config_number = 2^spin_number
 
+    adjacency_matrix = convert(Array{T,2}, gm)
+    prior_vector =  transpose(diag(adjacency_matrix))[1,:]
+
+    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, adj, prior, assignment_tmp) for i in (0:config_number-1)]
+    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)
@@ -51,14 +56,19 @@ function sample_generation{T <: Real}(samples_per_bin::Integer, adj::Array{T,2},
     return spin_samples
 end
 
-sample{T <: Real}(adjacency_matrix::Array{T,2}, number_sample::Integer) = sample(adjacency_matrix, number_sample, 1, Gibbs())[1]
-sample{T <: Real}(adjacency_matrix::Array{T,2}, number_sample::Integer, replicates::Integer) = sample(adjacency_matrix, number_sample, replicates, Gibbs())
+sample{T <: Real}(gm::FactorGraph{T}, number_sample::Integer) = sample(gm, number_sample, 1, Gibbs())[1]
+sample{T <: Real}(gm::FactorGraph{T}, number_sample::Integer, replicates::Integer) = sample(gm, number_sample, replicates, Gibbs())
 
-function sample{T <: Real}(adjacency_matrix::Array{T,2}, number_sample::Integer, replicates::Integer, sampler::Gibbs)
-    prior_vector = transpose(diag(adjacency_matrix))[1,:] #priors, or magnetic fields part
 
-    # generation of samples
-    samples = sample_generation(number_sample, adjacency_matrix, prior_vector, replicates)
+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)
 
     return samples
 end

test/common.jl

@@ -13,20 +13,20 @@ formulations = Dict(
 )
 
 gms = Dict(
-    "a" => [
+    "a" => FactorGraph([
         0.0 0.1 0.2;
         0.1 0.0 0.3;
         0.2 0.3 0.0
-    ],
-    "b" => [
+    ]),
+    "b" => FactorGraph([
         0.3 0.1 0.2;
         0.1 0.2 0.3;
         0.2 0.3 0.1
-    ],
-    "c" => [
+    ]),
+    "c" => FactorGraph([
         0.0 0.1 0.2 0.3;
         0.1 0.0 0.2 0.3;
         0.2 0.2 0.0 0.3;
         0.3 0.3 0.3 0.0
-    ]
+    ])
 )

test/data/build_data.jl

@@ -4,8 +4,8 @@
 
 using GraphicalModelLearning
 
-# remove old files
-`rm -rf *.csv`
+## remove old files
+#`rm -rf *.csv`
 
 include("../common.jl")
 

test/data/c_model_1.csv

@@ -0,0 +1,4 @@
+0.0, 0.1, 0.2, 0.3
+0.1, 0.0, 0.2, 0.3
+0.2, 0.2, 0.0, 0.3
+0.3, 0.3, 0.3, 0.0

test/data/c_model_2.csv

@@ -0,0 +1,4 @@
+0.1,  0.1, 0.2,  0.3
+0.0, -0.2, 0.2,  0.3
+0.0,  0.0, 0.3,  0.3
+0.0,  0.0, 0.0, -0.4

test/runtests.jl

@@ -5,6 +5,18 @@ using Base.Test
 include("common.jl")
 
 
+@testset "factor graphs" begin
+    for (name, gm) in gms
+        matrix = convert(Array{Float64,2}, gm)
+        gm2 = FactorGraph(matrix)
+        for key in keys(gm)
+            @test isapprox(gm[key], gm2[key])
+            @test isapprox(gm[key], matrix[key...])
+        end
+    end
+end
+
+
 @testset "gibbs sampler" begin
     for (name, gm) in gms
         srand(0) # fix random number generator
@@ -90,7 +102,7 @@ srand(0) # fix random number generator
                 sample_histo = sample(gm, act.samples)
                 #learned_gm = inverse_ising(sample_histo, method=act.formulation)
                 learned_gm = learn(sample_histo, act.formulation)
-                max_error = maximum(abs.(gm - learned_gm))
+                max_error = maximum(abs.(convert(Array{Float64,2}, gm) - learned_gm))
                 @test max_error <= act.threshold
             end
         end
@@ -101,10 +113,10 @@ end
 
 srand(0) # fix random number generator
 @testset "docs example" begin
-    model = [0.0 0.1 0.2; 0.1 0.0 0.3; 0.2 0.3 0.0]
+    model = FactorGraph([0.0 0.1 0.2; 0.1 0.0 0.3; 0.2 0.3 0.0])
     samples = sample(model, 100000)
     learned = learn(samples)
 
-    err = abs.(model - learned)
+    err = abs.(convert(Array{Float64,2}, model) - learned)
     @test maximum(err) <= 0.01
 end