Added an implementation of the Hungarian algorithm for solving the optimal

assignment problem (in the new max_cost_assignment() routine).

--HG--
extra : convert_revision : svn%3Afdd8eb12-d10e-0410-9acb-85c331704f74/trunk%404149

dlib/optimization.h

@@ -11,6 +11,7 @@
 #include "optimization/optimization_oca.h"
 #include "optimization/optimization_trust_region.h"
 #include "optimization/optimization_least_squares.h"
+#include "optimization/max_cost_assignment.h"
 
 #endif // DLIB_OPTIMIZATIOn_HEADER
 

dlib/optimization/max_cost_assignment.h

+// Copyright (C) 2011  Davis E. King (davis@dlib.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#ifndef DLIB_MAX_COST_ASSIgNMENT_H__
+#define DLIB_MAX_COST_ASSIgNMENT_H__
+
+#include "max_cost_assignment_abstract.h"
+#include "../matrix.h"
+#include <vector>
+#include <deque>
+
+namespace dlib
+{
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename EXP>
+    typename EXP::type assignment_cost (
+        const matrix_exp<EXP>& cost,
+        const std::vector<long>& assignment
+    )
+    {
+        DLIB_ASSERT(cost.nr() == cost.nc(),
+            "\t type assignment_cost(cost,assignment)"
+            << "\n\t cost.nr(): " << cost.nr()
+            << "\n\t cost.nc(): " << cost.nc()
+            << "\n\t min(assignment): " << min(vector_to_matrix(assignment)) 
+            << "\n\t max(assignment): " << max(vector_to_matrix(assignment)) 
+            );
+#ifdef ENABLE_ASSERTS
+        // can't call max on an empty vector. So put an if here to guard against it.
+        if (assignment.size() > 0)
+        {
+            DLIB_ASSERT(0 <= min(vector_to_matrix(assignment)) && max(vector_to_matrix(assignment)) < cost.nr(),
+                "\t type assignment_cost(cost,assignment)"
+                << "\n\t cost.nr(): " << cost.nr()
+                << "\n\t cost.nc(): " << cost.nc()
+                << "\n\t min(assignment): " << min(vector_to_matrix(assignment)) 
+                << "\n\t max(assignment): " << max(vector_to_matrix(assignment)) 
+                );
+        }
+#endif
+
+        typename EXP::type temp = 0;
+        for (unsigned long i = 0; i < assignment.size(); ++i)
+        {
+            temp += cost(i, assignment[i]);
+        }
+        return temp;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    namespace impl
+    {
+        template <typename U, long NR, long NC, typename MM, typename L>
+        inline void compute_slack(
+            const long x,
+            std::vector<U>& slack,
+            std::vector<long>& slackx,
+            const matrix<U,NR,NC,MM,L>& cost,
+            const std::vector<U>& lx,
+            const std::vector<U>& ly
+        )
+        {
+            for (long y = 0; y < cost.nc(); ++y)
+            {
+                if (lx[x] + ly[y] - cost(x,y) < slack[y])
+                {
+                    slack[y] = lx[x] + ly[y] - cost(x,y);
+                    slackx[y] = x;
+                }
+            }
+        }
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename U, long NR, long NC, typename MM, typename L>
+    std::vector<long> max_cost_assignment (
+        const matrix<U,NR,NC,MM,L>& cost
+    )                         
+    {
+        // This algorithm only works if the elements of the cost matrix can be reliably 
+        // compared using operator==. However, comparing for equality with floating point
+        // numbers is not a stable operation. So you need to use an integer cost matrix.
+        COMPILE_TIME_ASSERT(std::numeric_limits<U>::is_integer);
+        DLIB_ASSERT(cost.nr() == cost.nc(),
+            "\t std::vector<long> max_cost_assignment(cost)"
+            << "\n\t cost.nr(): " << cost.nr()
+            << "\n\t cost.nc(): " << cost.nc()
+            );
+
+        using namespace dlib::impl;
+        /*
+            I based the implementation of this algorithm on the description of the
+            Hungarian algorithm on the following websites:
+                http://www.math.uwo.ca/~mdawes/courses/344/kuhn-munkres.pdf
+                http://www.topcoder.com/tc?module=Static&d1=tutorials&d2=hungarianAlgorithm
+
+            Note that this is the fast O(n^3) version of the algorithm.
+        */
+
+        if (cost.size() == 0)
+            return std::vector<long>();
+
+        std::vector<U> lx, ly;
+        std::vector<long> xy;
+        std::vector<long> yx;
+        std::vector<char> S, T;
+        std::vector<U> slack;
+        std::vector<long> slackx;
+        std::vector<long> aug_path;
+
+
+
+
+        // Initially, nothing is matched. 
+        xy.assign(cost.nc(), -1);
+        yx.assign(cost.nc(), -1);
+        /*
+            We maintain the following invariant:
+                Vertex x is matched to vertex xy[x] and
+                vertex y is matched to vertex yx[y].
+
+                A value of -1 means a vertex isn't matched to anything.  Moreover,
+                x corresponds to rows of the cost matrix and y corresponds to the
+                columns of the cost matrix.  So we are matching X to Y.
+        */
+
+
+        // Create an initial feasible labeling.  Moreover, in the following
+        // code we will always have: 
+        //     for all valid x and y:  lx[x] + ly[y] >= cost(x,y)
+        lx.resize(cost.nc());
+        ly.assign(cost.nc(),0);
+        for (long x = 0; x < cost.nr(); ++x)
+            lx[x] = max(rowm(cost,x));
+
+        // Now grow the match set by picking edges from the equality subgraph until
+        // we have a complete matching.
+        for (long match_size = 0; match_size < cost.nc(); ++match_size)
+        {
+            std::deque<long> q;
+
+            // Empty out the S and T sets
+            S.assign(cost.nc(), false);
+            T.assign(cost.nc(), false);
+
+            // clear out old slack values
+            slack.assign(cost.nc(), std::numeric_limits<U>::max());
+            slackx.resize(cost.nc());
+            /*
+                slack and slackx are maintained such that we always
+                have the following (once they get initialized by compute_slack() below):
+                    - for all y:
+                        - let x == slackx[y]
+                        - slack[y] == lx[x] + ly[y] - cost(x,y)
+            */
+
+            aug_path.assign(cost.nc(), -1);
+
+            for (long x = 0; x < cost.nc(); ++x)
+            {
+                // If x is not matched to anything
+                if (xy[x] == -1)
+                {
+                    q.push_back(x);
+                    S[x] = true;
+
+                    compute_slack(x, slack, slackx, cost, lx, ly);
+                    break;
+                }
+            }
+
+
+            long x_start = 0;
+            long y_start = 0;
+
+            // Find an augmenting path.  
+            bool found_augmenting_path = false;
+            while (!found_augmenting_path)
+            {
+                while (q.size() > 0 && !found_augmenting_path)
+                {
+                    const long x = q.front();
+                    q.pop_front();
+                    for (long y = 0; y < cost.nc(); ++y)
+                    {
+                        if (cost(x,y) == lx[x] + ly[y] && !T[y])
+                        {
+                            // if vertex y isn't matched with anything
+                            if (yx[y] == -1) 
+                            {
+                                y_start = y;
+                                x_start = x;
+                                found_augmenting_path = true;
+                                break;
+                            }
+
+                            T[y] = true;
+                            q.push_back(yx[y]);
+
+                            aug_path[yx[y]] = x;
+                            S[yx[y]] = true;
+                            compute_slack(yx[y], slack, slackx, cost, lx, ly);
+                        }
+                    }
+                }
+
+                if (found_augmenting_path)
+                    break;
+
+
+                // Since we didn't find an augmenting path we need to improve the 
+                // feasible labeling stored in lx and ly.  We also need to keep the
+                // slack updated accordingly.
+                U delta = std::numeric_limits<U>::max();
+                for (unsigned long i = 0; i < T.size(); ++i)
+                {
+                    if (!T[i])
+                        delta = std::min(delta, slack[i]);
+                }
+                for (unsigned long i = 0; i < T.size(); ++i)
+                {
+                    if (S[i])
+                        lx[i] -= delta;
+
+                    if (T[i])
+                        ly[i] += delta;
+                    else
+                        slack[i] -= delta;
+                }
+
+
+
+                q.clear();
+                for (long y = 0; y < cost.nc(); ++y)
+                {
+                    if (!T[y] && slack[y] == 0)
+                    {
+                        // if vertex y isn't matched with anything
+                        if (yx[y] == -1)
+                        {
+                            x_start = slackx[y];
+                            y_start = y;
+                            found_augmenting_path = true;
+                            break;
+                        }
+                        else
+                        {
+                            T[y] = true;
+                            if (!S[yx[y]])
+                            {
+                                q.push_back(yx[y]);
+
+                                aug_path[yx[y]] = slackx[y];
+                                S[yx[y]] = true;
+                                compute_slack(yx[y], slack, slackx, cost, lx, ly);
+                            }
+                        }
+                    }
+                }
+            } // end while (!found_augmenting_path)
+
+            // Flip the edges along the augmenting path.  This means we will add one more
+            // item to our matching.
+            for (long cx = x_start, cy = y_start, ty; 
+                 cx != -1; 
+                 cx = aug_path[cx], cy = ty)
+            {
+                ty = xy[cx];
+                yx[cy] = cx;
+                xy[cx] = cy;
+            }
+
+        }
+
+
+        return xy;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+}
+
+#endif // DLIB_MAX_COST_ASSIgNMENT_H__
+
+

dlib/optimization/max_cost_assignment_abstract.h

+// Copyright (C) 2011  Davis E. King (davis@dlib.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#undef DLIB_MAX_COST_ASSIgNMENT_ABSTRACT_H__
+#ifdef DLIB_MAX_COST_ASSIgNMENT_ABSTRACT_H__
+
+#include "../matrix.h"
+#include <vector>
+
+namespace dlib
+{
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename EXP>
+    typename EXP::type assignment_cost (
+        const matrix_exp<EXP>& cost,
+        const std::vector<long>& assignment
+    );
+    /*!
+        requires
+            - cost.nr() == cost.nc()
+            - for all valid i:
+                - 0 <= assignment[i] < cost.nr()
+        ensures
+            - Interprets cost as a cost assignment matrix. That is, cost(i,j) 
+              represents the cost of assigning i to j.  
+            - Interprets assignment as a particular set of assignments. That is,
+              i is assigned to assignment[i].
+            - returns the cost of the given assignment. That is, returns
+              a number which is:
+                sum over i: cost(i,assignment[i])
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename U, long NR, long NC, typename MM, typename L>
+    std::vector<long> max_cost_assignment (
+        const matrix<U,NR,NC,MM,L>& cost
+    );
+    /*!
+        requires
+            - U == some integer type (e.g. int)
+            - cost.nr() == cost.nc()
+        ensures
+            - Finds and returns the solution to the following optimization problem:
+
+                Maximize: f(A) == assignment_cost(cost, A)
+                Subject to the following constraints:
+                    - The elements of A are unique. That is, there aren't any 
+                      elements of A which are equal.  
+                    - A.size() == cost.nr()
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+}
+
+#endif // DLIB_MAX_COST_ASSIgNMENT_ABSTRACT_H__
+

dlib/test/CMakeLists.txt

@@ -59,12 +59,13 @@ set (tests
    matrix_eig.cpp
    matrix_lu.cpp
    matrix_qr.cpp
+   max_cost_assignment.cpp
    md5.cpp
    member_function_pointer.cpp
    metaprogramming.cpp
    multithreaded_object.cpp
-   one_vs_one_trainer.cpp
    one_vs_all_trainer.cpp
+   one_vs_one_trainer.cpp
    optimization.cpp
    optimization_test_functions.cpp
    opt_qp_solver.cpp

dlib/test/makefile

@@ -69,6 +69,7 @@ SRC += matrix.cpp
 SRC += matrix_eig.cpp
 SRC += matrix_lu.cpp
 SRC += matrix_qr.cpp
+SRC += max_cost_assignment.cpp
 SRC += md5.cpp
 SRC += member_function_pointer.cpp
 SRC += metaprogramming.cpp

dlib/test/max_cost_assignment.cpp

+// Copyright (C) 2011  Davis E. King (davis@dlib.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+
+
+#include <dlib/optimization.h>
+#include <sstream>
+#include <string>
+#include <cstdlib>
+#include <ctime>
+#include <vector>
+#include "../rand.h"
+
+#include "tester.h"
+
+
+namespace  
+{
+
+    using namespace test;
+    using namespace dlib;
+    using namespace std;
+
+    logger dlog("test.max_cost_assignment");
+
+// ----------------------------------------------------------------------------------------
+
+    std::vector<std::vector<long> > permutations (
+        matrix<long,1,0> vals
+    )
+    {
+        if (vals.size() == 0)
+        {
+            return std::vector<std::vector<long> >();
+        }
+        else if (vals.size() == 1)
+        {
+            return std::vector<std::vector<long> >(1,std::vector<long>(1,vals(0)));
+        }
+
+
+        std::vector<std::vector<long> > temp;
+
+
+        for (long i = 0; i < vals.size(); ++i)
+        {
+            const std::vector<std::vector<long> >& res = permutations(remove_col(vals,i));       
+
+            for (unsigned long j = 0; j < res.size(); ++j)
+            {
+                temp.resize(temp.size()+1);
+                std::vector<long>& part = temp.back();
+                part.push_back(vals(i));
+                part.insert(part.end(), res[j].begin(), res[j].end());
+            }
+        }
+
+
+        return temp;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename T>
+    std::vector<long> brute_force_max_cost_assignment (
+        matrix<T> cost
+    )
+    {
+        if (cost.size() == 0)
+            return std::vector<long>();
+
+        const std::vector<std::vector<long> >& perms = permutations(range(0,cost.nc()-1));
+
+        T best_cost = std::numeric_limits<T>::min();
+        unsigned long best_idx = 0;
+        for (unsigned long i = 0; i < perms.size(); ++i)
+        {
+            const T temp = assignment_cost(cost, perms[i]);
+            if (temp > best_cost)
+            {
+                best_idx = i;
+                best_cost = temp;
+            }
+        }
+
+        return perms[best_idx];
+    }
+
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+
+    class test_max_cost_assignment : public tester
+    {
+    public:
+        test_max_cost_assignment (
+        ) :
+            tester ("test_max_cost_assignment",
+                    "Runs tests on the max_cost_assignment function.")
+        {}
+
+        dlib::rand::float_1a rnd;
+
+        template <typename T>
+        void test_hungarian()
+        {
+            long size = rnd.get_random_32bit_number()%7;
+            long range = rnd.get_random_32bit_number()%100;
+            matrix<T> cost = matrix_cast<T>(randm(size,size,rnd)*range) - range/2;
+
+            // use a uniform cost matrix sometimes
+            if ((rnd.get_random_32bit_number()%100) == 0)
+                cost = rnd.get_random_32bit_number()%100;
+
+            // negate the cost matrix every now and then
+            if ((rnd.get_random_32bit_number()%100) == 0)
+                cost = -cost;
+
+
+            std::vector<long> assign = brute_force_max_cost_assignment(cost);
+            const T true_eval = assignment_cost(cost, assign);
+            assign = max_cost_assignment(cost);
+
+            DLIB_TEST(assignment_cost(cost,assign) == true_eval);
+        }
+
+        void perform_test (
+        )
+        {
+            for (long i = 0; i < 1000; ++i)
+            {
+                if ((i%100) == 0)
+                    print_spinner();
+
+                test_hungarian<int>();
+                test_hungarian<long>();
+            }
+        }
+    } a;
+
+}
+
+
+