Structural Bioinformatics Library
Template C++ / Python API for developping structural bioinformatics applications.
User Manual

Earth_mover_distance

Authors: F. Cazals and T. Dreyfus and D. Mazauric

Introduction

Optimal transportation problems have a long standing history in mathematics and computer science, originating with the works of Monge on earth moving ( << la théorie des déblais et des remblais >> ) [88] . Such problems were later rephrased in terms of Riemannian geometry and measure theory [110] , one key concept being the distance between two distributions, namely the minimal amount of work that must be performed to transform one distribution into the other by moving distribution mass around.

Optimal transportation based metrics are now widely used, with applications in image processing [101] , statistics [77] , geometry processing [105] , biophysics [30] ...

In this context, this package provides two algorithms:

  • The classical i.e. linear program based Earth Mover Distance (EMD),

  • A version of EMD working on graphs rather than isolated sites, and accommodating connectivity constraints on these graphs.

Earth Mover Distance

We consider two sets of nodes, namely the source nodes $ \nodesS = \{ \nodesource{i} \}_{i=1,\dots, n_s} $, and the demand nodes $ \nodesD = \{ \nodedemand{j} \}_{j=1,\dots, n_d} $ respectively. We assume that these nodes are embedded in some metric space $ \calC $, and denote $ \dCalC(\cdot, \cdot) $ a metric in this space. We assume that each source node has a mass of supply (aka weights) denoted $ \sbasinw>0 $, and equivalently, that each demand node has a mass of demand $ \dbasinw $. The sums of masses for all source and demand nodes are denoted $ W_s $ and $ W_d $, that is $ \sum_{i =1,\dots,n_s} \sbasinw = W_s $ and $ \sum_{j =1,\dots,n_d} \dbasinw = W_d $. We also assume that transporting a unit mass of supply from $ \nodesource{i} $ to demand node $ \nodedemand{j} $ has unit cost $ \cij = \dCalC(\nodesource{i}, \nodedemand{j}) $.

A transport plan is a mathematical program aiming at satisfying the demand of the demand nodes given the supply of the supply nodes, while minimizing the overall transportation cost. More formally, a transport plan from $ \nodesS $ to $ \nodesD $ is defined by triples $ (\nodesource{i}, \nodedemand{j}, \flowij) $, with $ \flowij>0 $ the so-called flow quantities. Note that there are at most $ n_s\times n_d $ such triples.

Finding the optimal i.e. least cost transport plan amounts to solving the following linear program (LP):

\begin{equation} LP \begin{cases} \mbox{Min} \sum_{i=1,\dots,n_s, j=1,\dots,n_d} \flowij \times \cij\\ (1)\sum_{i =1,\dots,n_s} \flowij \leq \dbasinw & \forall j \in 1,\dots,n_d, \\ (2)\sum_{j = 1,\dots,n_d} \flowij \leq \sbasinw & \forall i \in 1,\dots,n_s, \\ (3)\flowij \geq 0 & \forall i \in 1,\dots,n_s, \forall j \in 1,\dots,n_d,\\ (4)\sum_{i =1,\dots,n_s} \sum_{j = 1,\dots,n_d} \flowij = \min(W_s, W_d). \end{cases} \end{equation}

The first equation is the linear functional to be minimized. Constraints (1,2) respectively state that a demand node does not receive more than its demand , and that a source node does not export more than it supply. Constraint (4) states that the total flow circulating is the minimum mass of supply or demand.

Based on this linear program, we introduce the total number of edges, the total flow, the total cost, and their ratio, known as the earth mover distance [101] :

\begin{equation} \begin{cases} \numedgesPemd = \sum_{i,j \mid \flowij > 0} 1,\\ \flowPemd = \sum_{i,j} \flowij,\\ \costPemd = \sum_{i,j} \flowij \cij,\\ \distPemd = \costPemd / \flowPemd. \end{cases} \end{equation}

Earth Mover Distance with connectivity constraints

The basic problem.

Consider now the case where the source and demand nodes are the vertices of two graphs $ G_S $ and $ G_D $. The solution provided by the EMD is oblivious to the connectivity of these two graphs. In selected applications, one may wish to preserve connectivity based properties. Typically, one would like a vertex or an edge of $ G_S $ to export flow to a connected component of $ G_D $. In the sequel, we specify a new problem such that such constraints are respected–which requires using more stringent transport plans.

More formally, to see whether a transport plan is valid, pick any connected subgraph $ H $ from $ \nodesS $. Let $ H' $ be the subgraph of $ G_D $ induced by the vertices from $ \nodesD $ such that for each such vertex $ d_j $, there is at least one edge emanating from a vertex $ s_i $ of the subgraph $ H $ with $ \flowij>0 $. The transport plan is called valid iff the subgraph $ I_D $ is also connected.

We summarize this discussion with the following definition:

A transport plan is said to respect connectivity constraints iff any connected set of vertices from $ \nodesS $ induces, through edges carrying strictly positive flow, a connected set of vertices from $ \nodesD $.


The following remarks are in order:

  • A transport plan respecting connectivity constraints may not fully satisfy the demand.

  • It is actually easy to respect all connectivity constraints, by using an arbitrary number of edges which may carry an arbitrarily small amount of flow, and guarantee that connectivity constraints are respected [26] . This motivates the definition of a new problem denoted $ \text{\pbemdccc} $.

The problem solved: $ \text{\pbemdccc} $.

An alternative to the ill-posed problem of minimizing the transport cost while conserving connectivity constraints is the so-called maximum-flow under cost and connectivity constraints problem ( $ \text{\pbemdccc} $). In this problem one aims at computing the largest amount of flow that can be supported respecting the connectivity constraints and such that the total cost is less than a given bound $ C $. We define the following upper bound for the maximum total cost of any admissible flow as

\begin{equation} C_{max} = \sum_{j=1,\dots,n_d} \dbasinw \max_{i=1,\dots,n_s} \cij. \end{equation}

Consider an upper bound $ C \leq C_{max} $, and an upper bound $ 0 \leq M \leq n_s\times n_d $ on the number of edges which carry non null flow. Denote $ \mathcal{H} $ and $ \mathcal{H'} $ the set of all connected subgraphs of $ G_S $ and $ G_D $, respectively. Problem $ \text{\pbemdccc} $ is defined as follows:

\begin{equation} \left\{ \begin{array}{ll} \mbox{Maximize} \sum_{i=1,\dots,n_s} \sum_{j=1,\dots,n_d} \flowij \\ \mbox{subject to} \\ (1)\sum_{i=1,\dots,n_s} \flowij \leq \dbasinw ~~~ \forall j=1,\dots,n_d, \\ (2)\sum_{j=1,\dots,n_d} \flowij \leq \sbasinw ~~~ \forall i=1,\dots,n_s, \\ (3)\flowij \geq 0 ~~~ \forall i=1,\dots,n_s, \forall j=1,\dots,n_d, \\ (4)\sum_{i=1,\dots,n_s} \sum_{j=1,\dots,n_d} \flowij ~ c_{i,j} \leq C, \\ (5)\sum_{i=1,\dots,n_s} \sum_{j=1,\dots,n_d \mid \flowij > 0} 1 \leq M, \\ (6)H' = \{v'_{j} \mid \flowij > 0, v_{i} \in H\} \in \mathcal{H'} ~~~ \forall H \in \mathcal{H} \end{array} \right. \end{equation}

Constraints (1,2,3) are the usual constraints found in EMD. Constraint (4) takes into account the upper bound on the transport cost, while constraint (5) bounds the number of edges used. Finally, constraint (6) encodes the connectivity constraints.

Algorithms

Earth Mover Distance

Algorithms.

Once the weights of all vertices, solving the linear program of Eq. (emd-eq-emd-lp) has polynomial complexity [69] . Practically, various solvers can be used, e.g. the one from the Computational Geometry Algorithms Library [32] , lp_solve, the CPLEX solver from IBM, etc. The algorithm used in this package is lp_solve, which take as input the MPS file format.

In the following, the algorithm solving the linear program of Eq. (emd-eq-emd-lp) is called $ \text{\algoEMDLP} $.

Earth Mover Distance with connectivity constraints

Algorithms.

Finding transport plans respecting connectivity constraints are hard combinatorial problems. The following results are proved in [26] :

  • the decision version of $ \text{\pbemdccc} $ is NP-complete,

  • a Polynomial Time Approximation Scheme for $ \text{\pbemdccc} $ exists when transport plans involving a large number of edges are allowed.

Practically, we solve $ \text{\pbemdccc} $ with two algorithms denoted $ \text{\algoEMDCCCG} $ and $ \text{\algoEMDCCCGI} $. These algorithms are implemented within the class SBL::CADS::T_Earth_mover_distance_with_connectivity_constraints , and can be called by specifying a recursion mode:

  • Recursion mode norec: algorithm $ \text{\algoEMDCCCG} $. A greedy strategy returning an admissible solution, given a maximum cost $ C_{max} $ provided by the user.

  • Recursion mode refined: algorithm $ \text{\algoEMDCCCGI} $. Given a range $ [0, C_{max}] $ of costs, this interval is iteratively halved, and the upper bound of each interval is used to compute a solution with this upper bound as constraint.

  • Recursion mode coarse: algorithm $ \text{\algoEMDCCCGI} $. As in the previous option, except that upon splitting the interval into two, only the left interval (that associated with the smallest upper bound) is used to recurse.

For a given solution, in a manner analogous to Eq. (emd-eq-emd-lp-dist), we define the total number of edges, the total flow, the total cost, and their ratio:

\begin{equation} \begin{cases} \numedgesAemdccc = \sum_{i,j \mid \flowij > 0} 1,\\ \flowAemdccc = \sum_{i,j} \flowij,\\ \costAemdccc = \sum_{i,j} \flowij \cij,\\ \distAemdccc = \costAemdccc / \flowAemdccc. \end{cases} \end{equation}

It can be shown that the Earth Mover Distance as defined by Eq. (emd-eq-emd-lp) yields a metric, provided that $ \dCalC $ is itself a metric, and that the sum of weights for the source and the demand graphs are equal [101] . On the other hand, earth mover distances with connectivity constraints are not symmetric in general [26] , but can easily be made so by taking a symmetric function of the one sided quantities. To assess this lack of symmetry, we introduce the following ratios, respectively geared towards the flow and the cost:

\begin{equation} \begin{cases} \ratiosymFlow = \frac{ \min(\flowAemdccc(A, B), \flowAemdccc(B, A)) } { \max(\flowAemdccc(A, B), \flowAemdccc(B, A))} ,\\ \\ \ratiosymCost = \frac{ \min(\costAemdccc(A, B), \costAemdccc(B, A)) } {\max(\costAemdccc(A, B), \costAemdccc(B, A))}. \end{cases} \end{equation}


Implementation and functionality

Vertices information.

The generic paradigm followed in the SBL implies that any kind of data structure could be used as input of the EMD algorithms. In particular, it could be two ensembles of geometric points, or two graphs where vertices are embedded in a metric space. For this reason, the algorithms are templated by a concept VerticesAccessor defining how to access the data used in the algorithms. (NB: an accessor allows iterating over the vertices and their neighbors, if any, and also allows access to the mass of vertices.) Two models of VerticesAccessor are predefined:

EMD algorithms.

Each algorithm is implemented in a separate class as a functor taking as argument the source and the demand, and returning the output transportation plan:

  • SBL::CADS::T_Earth_mover_distance < VerticesAccessor , Distance > is the EMD algorithm described in section Earth Mover Distance . The template parameter VerticesAccessor is described just above, and Distance is a distance functor between the objects represented by the vertices.

Connectivity checker.

In order to compare both algorithms described in sections Earth Mover Distance and Earth Mover Distance with connectivity constraints , the class SBL::CADS::T_Earth_mover_distance_connectivity_constraints_checker< EarthMoverDistance > computes the number of edges in a transportation plan that do not respect the connectivity of the demand. The template parameter EarthMoverDistance is one of the EMD algorithms class. Note that for a transportation plan computed from the class SBL::CADS::T_Earth_mover_distance_with_connectivity_constraints, the connectivity must be respected.

Module.

This package provides a module, see Modules, to be integrated into a workflow : SBL::CADS::T_Earth_mover_distance_module< ModuleTraits > . The template parameter ModuleTraits should define two types:

  • EMD_vertices_accessor that has to comply the concept VerticesAccessor
  • EMD_distance that has to comply with the concept Distance .

This module allows adding options to select the algorithm used, to run symmetric instances of the problems (source and demand are reversed), etc.

In addition to the log file containing statistics on the cost and the flow, two output files are reported :

  • The transportation plan, as a Graphviz file, suffixed by __transportation_plan.dot : this file can be processed by Graphviz to produce a visual representation of the transportation plan as a bipartite graph: the vertices shown are the source and demand nodes; for each edge, the flow and cost are represented.

  • A serialized XML archive suffixed by __emd_engine.xml : the archive contains in particular the serialization of the transportation plan, to be used with the package PALSE for further analysis. In particular, for each edge of the transportation plan, the flow and associated cost are provided.

Examples

Earth Mover Distance

The following example computes the transportation plan of two sets of three weighted 2D points with no connectivity, using the algorithm described in section Earth Mover Distance , and then dumps statistics about the transportation plan.

#include <SBL/CADS/Earth_mover_distance.hpp>
#include <SBL/CADS/Earth_mover_distance_vertices_accessor_vector.hpp>
#include <CGAL/Cartesian.h>
typedef CGAL::Cartesian<double>::Point_2 Point_2;
<Point_2, double> EMD_vertices_accessor;
class Distance
{
public:
typedef double FT;
inline FT operator()(const Point_2& u, const Point_2& v)const
{
return CGAL::sqrt(CGAL::squared_distance(u, v));
}
};
int main(int argc, char *argv[])
{
EMD::set_verbose_mode(1);
std::vector<EMD_vertices_accessor::Weighted_vertex> source;
source.push_back(EMD_vertices_accessor::Weighted_vertex(Point_2(0, 0), 4));
source.push_back(EMD_vertices_accessor::Weighted_vertex(Point_2(1, 0), 2));
source.push_back(EMD_vertices_accessor::Weighted_vertex(Point_2(2, 0), 4));
std::vector<EMD_vertices_accessor::Weighted_vertex> demand;
demand.push_back(EMD_vertices_accessor::Weighted_vertex(Point_2(0, 0), 2.5));
demand.push_back(EMD_vertices_accessor::Weighted_vertex(Point_2(0, 1), 5));
demand.push_back(EMD_vertices_accessor::Weighted_vertex(Point_2(0, 2), 2.5));
EMD emd;
EMD::Transportation_plan tplan = emd(source, demand);
std::cout << std::endl << tplan << std::endl << std::endl;
std::cout << "Source flow: " << tplan.get_final_flow() << " / " << tplan.get_total_source() << " = " << tplan.get_percentage_of_used_flow() << std::endl;
std::cout << "Remaining demand: " << tplan.get_remaining_demand() << " / " << tplan.get_total_demand() << " = " << tplan.get_percentage_of_remaining_demand() << std::endl;
std::cout << "Ratio flow/demand: " << tplan.get_demand_flow_ratio() << std::endl;
std::cout << "Transportation cost (normalized): " << tplan.get_transportation_cost() << " / " << tplan.get_final_flow() << " = " << tplan.get_normalized_transportation_cost() << std::endl;
return 0;
}

Earth Mover Distance with connectivity constraints

The following example computes the transportation plan of two sets of three weighted 2D points with no connectivity, using the algorithm described in section Earth Mover Distance with connectivity constraints , and then dumps statistics about the transportation plan.

#include <SBL/CADS/Earth_mover_distance_with_connectivity_constraints.hpp>
#include <SBL/CADS/Earth_mover_distance_vertices_accessor_graph.hpp>
#include <CGAL/Cartesian.h>
typedef CGAL::Cartesian<double>::Point_2 Point_2;
<Point_2, double> EMD_vertices_accessor;
class Distance
{
public:
typedef double FT;
inline FT operator()(const Point_2& u, const Point_2& v)const
{
return CGAL::sqrt(CGAL::squared_distance(u, v));
}
};
typedef EMD_vertices_accessor::Vertices_container Graph;
int main(int argc, char *argv[])
{
EMD::set_verbose_mode(1);
Graph source;
boost::add_vertex(EMD_vertices_accessor::Weighted_vertex(Point_2(0, 0), 4), source);
boost::add_vertex(EMD_vertices_accessor::Weighted_vertex(Point_2(1, 0), 2), source);
boost::add_vertex(EMD_vertices_accessor::Weighted_vertex(Point_2(2, 0), 4), source);
Graph demand;
boost::add_vertex(EMD_vertices_accessor::Weighted_vertex(Point_2(0, 0), 2.5), demand);
boost::add_vertex(EMD_vertices_accessor::Weighted_vertex(Point_2(0, 1), 5), demand);
boost::add_vertex(EMD_vertices_accessor::Weighted_vertex(Point_2(0, 2), 2.5), demand);
std::cout << "Modes: " << EMD::get_algorithm_type() << ", " << EMD::get_edge_selection_strategy() << ", " << EMD::get_recursion_mode() << std::endl;
EMD emd;
EMD::Transportation_plan tplan = emd(source, demand);
std::cout << std::endl << tplan << std::endl << std::endl;
std::cout << "Source flow: " << tplan.get_final_flow() << " / " << tplan.get_total_source() << " = " << tplan.get_percentage_of_used_flow() << std::endl;
std::cout << "Remaining demand: " << tplan.get_remaining_demand() << " / " << tplan.get_total_demand() << " = " << tplan.get_percentage_of_remaining_demand() << std::endl;
std::cout << "Ratio flow/demand: " << tplan.get_demand_flow_ratio() << std::endl;
std::cout << "Transportation cost (normalized): " << tplan.get_transportation_cost() << " / " << tplan.get_final_flow() << " = " << tplan.get_normalized_transportation_cost() << std::endl;
return 0;
}

Earth Mover Distance connectivity constraints checker

The following example does the same calculations as in section Earth Mover Distance, but also checks that the connectivity constraints are respected. If not, the number of edges that violate the constraints is dumped.

#include <SBL/CADS/Earth_mover_distance.hpp>
#include <SBL/CADS/Earth_mover_distance_connectivity_constraints_checker.hpp>
#include <SBL/CADS/Earth_mover_distance_vertices_accessor_graph.hpp>
#include <CGAL/Cartesian.h>
typedef CGAL::Cartesian<double>::Point_2 Point_2;
<Point_2, double> EMD_vertices_accessor;
class Distance
{
public:
typedef double FT;
inline FT operator()(const Point_2& u, const Point_2& v)const
{
return CGAL::sqrt(CGAL::squared_distance(u, v));
}
};
typedef EMD_vertices_accessor::Vertices_container Graph;
int main(int argc, char *argv[])
{
Graph source;
boost::add_vertex(EMD_vertices_accessor::Weighted_vertex(Point_2(0, 0), .4), source);
boost::add_vertex(EMD_vertices_accessor::Weighted_vertex(Point_2(1, 0), .2), source);
boost::add_vertex(EMD_vertices_accessor::Weighted_vertex(Point_2(2, 0), .4), source);
Graph demand;
boost::add_vertex(EMD_vertices_accessor::Weighted_vertex(Point_2(0, 0), .25), demand);
boost::add_vertex(EMD_vertices_accessor::Weighted_vertex(Point_2(0, 1), .5), demand);
boost::add_vertex(EMD_vertices_accessor::Weighted_vertex(Point_2(0, 2), .25), demand);
EMD::set_verbose_mode(1);
EMD emd;
EMD::Transportation_plan tplan = emd(source, demand);
EMD_cc_checker::set_verbose_mode(1);
EMD_cc_checker checker;
EMD_cc_checker::result_type res = checker(source, demand, tplan);
tplan.dump(std::cout);
std::cout << "vertices fraction : " << res.first << std::endl;
std::cout << "edges fraction : " << res.second << std::endl;
return 0;
}

Earth Mover Distance module

The following example shows a simple instantiation of the EMD module within a workflow.

#include <SBL/CADS/Earth_mover_distance_vertices_accessor_graph.hpp>
#include <SBL/Modules/Earth_mover_distance_module.hpp>
#include <SBL/Modules/Module_based_workflow.hpp>
#include <CGAL/Cartesian.h>
typedef CGAL::Cartesian<double>::Point_2 Point_2;
class Traits
{
public:
<Point_2, double> EMD_vertices_accessor;
class EMD_distance
{
public:
typedef double FT;
inline FT operator()(const Point_2& u, const Point_2& v)const
{
return CGAL::sqrt(CGAL::squared_distance(u, v));
}
};//end class EMD_distance
};//end class Traits
typedef Traits::EMD_vertices_accessor::Vertices_container Graph;
class Workflow
{
public:
private:
private:
Workflow_vertex m_EMD;
EMD_module* m_EMD_module;
public:
inline Workflow()
: m_workflow("example_emd_module")
{
this->m_workflow.add_helper("Example of instantiation of the EMD module.");
this->m_EMD = this->m_workflow.register_module<EMD_module>("EMD");
this->m_EMD_module = (EMD_module*)&this->m_workflow.get_module(this->m_EMD);
this->m_workflow.set_start_module(this->m_EMD);
this->m_workflow.set_end_module(this->m_EMD);
}
inline void start(int argc, char *argv[])
{
this->m_workflow.parse_command_line(argc, argv);
this->m_EMD_module->get_source() = new Graph();
boost::add_vertex(Traits::EMD_vertices_accessor::Weighted_vertex(Point_2(0, 0), .4),
*this->m_EMD_module->get_source());
boost::add_vertex(Traits::EMD_vertices_accessor::Weighted_vertex(Point_2(1, 0), .2),
*this->m_EMD_module->get_source());
boost::add_vertex(Traits::EMD_vertices_accessor::Weighted_vertex(Point_2(2, 0), .4),
*this->m_EMD_module->get_source());
this->m_EMD_module->get_demand() = new Graph();
boost::add_vertex(Traits::EMD_vertices_accessor::Weighted_vertex(Point_2(0, 0), .25),
*this->m_EMD_module->get_demand());
boost::add_vertex(Traits::EMD_vertices_accessor::Weighted_vertex(Point_2(0, 1), .5),
*this->m_EMD_module->get_demand());
boost::add_vertex(Traits::EMD_vertices_accessor::Weighted_vertex(Point_2(0, 2), .25),
*this->m_EMD_module->get_demand());
this->m_workflow.start();
}
};//end class Workflow
int main(int argc, char *argv[])
{
Workflow workflow;
workflow.start(argc, argv);
return 0;
}

Applications

Practically, the sources of the examples above are provided in the examples section.

In the context of biophysics the package Energy_landscape_comparison uses the Earth_mover_distance algorithms to compare transition graphs (see Transition_graph_traits) derived from sampled energy landscapes.

In addition, this package defines four programs for using EMD and EMD CC in different contexts:

  • sbl-emd-wp-euclid.exe : a program comparing two ensembles of weighted D-dimensional points, using the Euclidean metric; the input consists of two text files listing the points (dimension followed by coordinates, one point per line), and two text files listing the weights of the points. (NB: the ordering of points and weight must be the same.)
  • sbl-emd-graph-euclid.exe : a program comparing vertex-weighted graphs where vertices are D-dimensional points, using the Euclidean metric; the input consists on two text files listing the points (dimension followed by coordinates, one point per line), two text files listing the weights of the points in the order they are listed, and two text files listing the edges of the graph (positions of points in their file, starting at 0, one edge per line). For example, this program can be used to compare the critical points of two elevated surfaces as defined in the package Morse_theory_based_analyzer : the program sbl-Morse-theory-based-analyzer-nng-euclid.exe used with the option –report-MSW-as-txt reports exactly the input of sbl-emd-graph-euclid.exe.
  • sbl-emd-tg-lrmsd.exe : a program comparing two transition graphs using the lRMSD metric, as defined in the package Energy_landscape_comparison.