Added the find_max_factor_graph_viterbi() routine for performing MAP

inference in chain-structured factor graphs.

dlib/optimization.h

@@ -14,6 +14,7 @@
 #include "optimization/max_cost_assignment.h"
 #include "optimization/max_sum_submatrix.h"
 #include "optimization/find_max_factor_graph_nmplp.h"
+#include "optimization/find_max_factor_graph_viterbi.h"
 
 #endif // DLIB_OPTIMIZATIOn_HEADER
 

dlib/optimization/find_max_factor_graph_viterbi.h

+// Copyright (C) 2011  Davis E. King (davis@dlib.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#ifndef DLIB_FIND_MAX_FACTOR_GRAPH_VITERBi_H__
+#define DLIB_FIND_MAX_FACTOR_GRAPH_VITERBi_H__
+
+#include "find_max_factor_graph_viterbi_abstract.h"
+#include <vector>
+#include "../matrix.h"
+#include "../array2d.h"
+
+
+namespace dlib
+{
+
+// ----------------------------------------------------------------------------------------
+
+    namespace impl
+    {
+        struct viterbi_data
+        {
+            viterbi_data() :val(0), back_index(0) {}
+            double val;
+            unsigned long back_index;
+        };
+
+        template <long NC>
+        inline bool advance_state(
+            matrix<unsigned long,1,NC>& node_states,
+            unsigned long num_states
+        )
+        /*!
+            ensures
+                - advances node_states to the next state by adding 1
+                  to node_states(node_states.size()-1) and carrying any 
+                  rollover (modulo num_states).  Stores the result into #node_states.
+                - if (#node_states is all zeros) then
+                    - returns false
+                - else
+                    - returns true
+        !*/
+        {
+            for (long i = node_states.size()-1; i >= 0; --i)
+            {
+                node_states(i) += 1;
+                if (node_states(i) < num_states)
+                    return true;
+
+                node_states(i) = 0;
+            }
+            return false;
+        }
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename map_problem
+        >
+    void find_max_factor_graph_viterbi (
+        const map_problem& prob,
+        std::vector<unsigned long>& map_assignment
+    )
+    {
+        using namespace dlib::impl;
+        const unsigned long order = map_problem::order;
+        const unsigned long num_states = map_problem::num_states;
+        COMPILE_TIME_ASSERT(num_states > 0);
+
+
+        DLIB_ASSERT(std::pow(num_states,order) < std::numeric_limits<unsigned long>::max(),
+            "\t void find_max_factor_graph_viterbi()"
+            << "\n\t The order is way too large for this algorithm to handle."
+            << "\n\t order:      " << order
+            << "\n\t num_states: " << num_states 
+            << "\n\t std::pow(num_states,order):                " << std::pow(num_states,order) 
+            << "\n\t std::numeric_limits<unsigned long>::max(): " << std::numeric_limits<unsigned long>::max() 
+            );
+
+        const unsigned long trellis_size = static_cast<unsigned long>(std::pow(num_states,order)); 
+        unsigned long init_ring_size = 1; 
+
+        array2d<impl::viterbi_data> trellis;
+        trellis.set_size(prob.number_of_nodes(), trellis_size);
+
+
+        for (unsigned long node = 0; node < prob.number_of_nodes(); ++node)
+        {
+
+            if (node < order)
+            {
+                matrix<unsigned long,1,0> node_states;
+                node_states.set_size(std::min<int>(node, order) + 1);
+                node_states = 0;
+
+                unsigned long idx = 0;
+                if (node == 0)
+                {
+                    do 
+                    {
+                        trellis[node][idx].val = prob.factor_value(node,node_states);
+                        ++idx;
+                    } while(advance_state(node_states,num_states));
+                }
+                else
+                {
+                    init_ring_size *= num_states;
+                    do 
+                    {
+                        const unsigned long back_index = idx%init_ring_size;
+                        trellis[node][idx].val = prob.factor_value(node,node_states) + trellis[node-1][back_index].val;
+                        trellis[node][idx].back_index = back_index;
+                        ++idx;
+                    } while(advance_state(node_states,num_states));
+
+                }
+            }
+            else
+            {
+                matrix<unsigned long,1,order+1> node_states;
+                node_states = 0;
+
+                unsigned long count = 0;
+                for (unsigned long i = 0; i < trellis_size; ++i)
+                {
+                    unsigned long back_index = 0;
+                    double best_score = -std::numeric_limits<double>::infinity();
+                    for (unsigned long s = 0; s < num_states; ++s)
+                    {
+                        const double temp = prob.factor_value(node,node_states) + trellis[node-1][count%trellis_size].val;
+                        if (temp > best_score)
+                        {
+                            best_score = temp;
+                            back_index = count%trellis_size;
+                        }
+
+                        advance_state(node_states,num_states);
+                        ++count;
+                    }
+
+                    trellis[node][i].val = best_score;
+                    trellis[node][i].back_index = back_index;
+                }
+            }
+        }
+
+
+        map_assignment.resize(prob.number_of_nodes());
+        // Figure out which state of the last node has the biggest value. 
+        unsigned long back_index = 0;
+        double best_val = -std::numeric_limits<double>::infinity();
+        for (long i = 0; i < trellis.nc(); ++i)
+        {
+            if (trellis[trellis.nr()-1][i].val > best_val)
+            {
+                best_val = trellis[trellis.nr()-1][i].val;
+                back_index = i;
+            }
+        }
+        // Follow the back links to find the decoding.
+        for (long node = map_assignment.size()-1; node >= 0; --node)
+        {
+            map_assignment[node] = back_index/init_ring_size;
+            back_index = trellis[node][back_index].back_index;
+            if (node < (long)order)
+                init_ring_size /= num_states;
+        }
+
+    }
+
+// ----------------------------------------------------------------------------------------
+
+}
+
+#endif // DLIB_FIND_MAX_FACTOR_GRAPH_VITERBi_H__
+
+

dlib/optimization/find_max_factor_graph_viterbi_abstract.h

+// Copyright (C) 2011  Davis E. King (davis@dlib.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#undef DLIB_FIND_MAX_FACTOR_GRAPH_VITERBi_ABSTRACT_H__
+#ifdef DLIB_FIND_MAX_FACTOR_GRAPH_VITERBi_ABSTRACT_H__
+
+#include <vector>
+#include "../matrix.h"
+
+
+namespace dlib
+{
+
+// ----------------------------------------------------------------------------------------
+
+    class map_problem 
+    {
+        /*!
+            WHAT THIS OBJECT REPRESENTS
+                This object represents a chain-structured factor graph or graphical 
+                model.  In particular, this object defines the interface a MAP problem 
+                on a factor graph must implement if it is to be solved using the 
+                find_max_factor_graph_viterbi() routine defined at the bottom of this file.  
+
+                Note that there is no dlib::map_problem object.  What you are looking 
+                at here is simply the interface definition for a map problem.  You must 
+                implement your own version of this object for the problem you wish to 
+                solve and then pass it to the find_max_factor_graph_viterbi() routine.
+        !*/
+
+    public:
+
+        // This model can represent a high order Markov chain.  If order==1 then map_problem
+        // represents a basic chain-structured graph where nodes only depend on their immediate
+        // neighbors.  However, high order Markov models can also be used by setting order > 1.
+        const static unsigned long order;
+
+        // Defines the number of states attainable by each variable/node in the graph.
+        const static unsigned long num_states;  
+
+        unsigned long number_of_nodes (
+        ) const;
+        /*!
+            ensures
+                - returns the number of nodes in the factor graph.  Or in other words, 
+                  returns the number of variables in the MAP problem.
+        !*/
+
+        template <
+            typename EXP 
+            >
+        double factor_value (
+            unsigned long node_id,
+            const matrix_exp<EXP>& node_states
+        ) const;
+        /*!
+            requires
+                - EXP::type == unsigned long
+                  (i.e. node_states contains unsigned longs)
+                - node_id < number_of_nodes()
+                - node_states.size() == min(node_id, order) + 1
+                - max(node_states) < num_states
+            ensures
+                - In a chain-structured graph, each node has a potential function associated with
+                  it.  The potential function operates on the variable given by the node as well
+                  as the order previous variables.  Therefore, factor_value() returns the value 
+                  of the factor/potential function associated with node node_id where the following 
+                  nodes take on the values defined below:
+                    - node_states(0) == the value of the node with ID node_id
+                    - node_states(i) == the value of the node with ID node_id-i
+        !*/
+
+    };
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename map_problem
+        >
+    void find_max_factor_graph_viterbi (
+        const map_problem& prob,
+        std::vector<unsigned long>& map_assignment
+    );
+    /*!
+        requires
+            - map_problem::num_states > 0
+            - std::pow(map_problem::num_states, map_problem::order) < std::numeric_limits<unsigned long>::max()
+              (i.e. The Viterbi algorithm is exponential in the order of the map problem.  So don't 
+              make order too large.)
+            - map_problem == an object with an interface compatible with the map_problem
+              object defined at the top of this file.
+        ensures
+            - This function is a tool for exactly solving the MAP problem in a chain-structured 
+              graphical model or factor graph.  That is, it attempts to solve a certain kind of 
+              optimization problem which can be defined as follows:
+                - Let X be a set of prob.number_of_nodes() integer valued variables, each taking
+                  a value in the range [0, map_problem::num_states).
+                - Let X(i) = the ith variable in X.
+                - Let F(i) = factor_value_i(X(i), X(i-1), ..., X(i-map_problem::order))
+                  (This is the value returned by prob.factor_value(i, node_states).  Note that
+                  each factor value function operates on at most map_problem::order+1 variables.
+                  Moreover, the variables are adjacent and hence the graph is "chain-structured".)
+
+                 Then this function finds the assignments to the X variables which  
+                 maximizes: sum over all valid i: F(i)
+
+            - #map_assignment == the result of the optimization.   
+            - #map_assignment.size() == prob.number_of_nodes()
+            - for all valid i:
+                - #map_assignment[i] < map_problem::num_states
+                - #map_assignment[i] == The MAP assignment for node/variable i.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+}
+
+#endif // DLIB_FIND_MAX_FACTOR_GRAPH_VITERBi_ABSTRACT_H__
+
+
+