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 >> ) [109] . Such problems were later rephrased in terms of Riemannian geometry and measure theory [138] , 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 [125] , statistics [95] , geometry processing [131] , biophysics [34] ...

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 [125] :

\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 [28] . 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 [84] . Practically, various solvers can be used, e.g. the one from the Computational Geometry Algorithms Library [37] , 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 [28] :

  • 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 [125] . On the other hand, earth mover distances with connectivity constraints are not symmetric in general [28] , 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

Programs

Executables in this package are named sbl-emd-OBJ-COST.exe, with

  • OJB : weighted point, graph: graph, tg: transition graph.
  • COST : the cost function used between points or vertices of the graphs compared.

Objects. The rationale for handling these three types of objects is as follows:

  • Weighted points: a weighted point is a pair (point in a metric space, weight). See the example below. Note that several metrics/costs may be used for the same point type.
  • Graph : a simple graph specified in text format. That is, assuming vertices are indexed in the range 0..n-1, edges are pairs of integers listed in a text file.
  • Transition graph: a particular graph type whose edges correspond to saddle points on (potential) energy landscapes, see the pacakge Transition_graph_traits

Costs. A function assigning a positive value to two points/graph vertices. Recall that is this function is a metric (satisfies the triangle inequality), then EMD is also a metric.

Examples.

Example 1: Consider d-dimensional points in the Euclidean space, for which one uses the Euclidean distance. We provide the three programs, detailed below:

  • sbl-emd-wp-euclid.exe, sbl-emd-graph-euclid.exe, sbl-emd-tg-euclid.exe

Example 2: Assume now that the d-dimensional points represent molecular conformations, for which one uses the least Root Mean Square Deviation as a distance. This package provides

Example 3: Assume now that the objects compared are sequences of nucleotides (nucleic acids) or amino-acids (polypeptide chains).

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.


The program sbl-emd-graph-euclid.exe 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.


The program sbl-emd-tg-euclid.exe compares two transition graphs, using the Euclidean metric, as defined in the package Energy_landscape_comparison. For example, the application Transition_graph_of_energy_landscape_builders creates such transition graphs from samples of the energy landscape of target molecule.


The program sbl-emd-tg-lrmsd.exe compares two transition graphs using the lRMSD metric, as defined in the package Energy_landscape_comparison.


Main specifications

Main options – comparing weighted points.

points1: Text file listing the first D-dimensional points weights1: Text file listing the first weights

The main options are:

–points-file string: Input points as D-dimensional points – passed twice
–weights-file string: Input weights – passed twice


Main options – comparing graphs or transition graphs. The following options, each called twice i.e. once for each graph, are used to passe the input:

–points-file string: Description of points in point d format – one point per line
–weights-file string: Description of positive weights – the ordering as points
–edges-file string: File listing edges of the graph – two integers in 0..n-1


Main options for EMD with connectivity constraints. The main options are:

–constraints string: use EMD-CC
–algorithm string: with EMD-CC –select the algorithm (reg or star)
–edgeSelection string: with EMD-CC – selection scheme (min-cost or first-found)
–recursion string: with EMD-CC – recursion mode (norec / refined / coarse)


Jupyter demo

See the following jupyter notebook:

  • Jupyter notebook file
  • Earth_mover_distance

    Earth_mover_distance

    Example : Comparing weighted points

    Options.

    The main options of the compareSamples method in the next cell are:

    • points1: Text file listing the first D-dimensional points
    • weights1: Text file listing the first weights
    • points2: Text file listing the second D-dimensional points
    • weights2: Text file listing the second weights
    In [1]:
    import re  #regular expressions
    import sys #misc system
    import os
    import pdb
    import shutil # python 3 only
    
    def compareSamples(points1, weights1, points2, weights2): 
        
        odir = "tmp-results-samples"
        if os.path.exists(odir):
            os.system("rm -rf %s" % odir)
        os.system( ("mkdir %s" % odir) )
        
        # check executable exists and is visible
        exe = shutil.which("sbl-emd-wp-euclid.exe")
        if exe:
            print(("Using executable %s\n" % exe))
            cmd = "sbl-emd-wp-euclid.exe --points-file %s --points-file %s -v -l -d %s" %\
                  (points1, points2, odir)
            if weights1 != None and weights2 != None:
                cmd += " --weights-file %s  --weights-file %s" % (weights1, weights2)
            os.system(cmd)
            
            cmd = "ls %s" % odir
            ofnames = os.popen(cmd).readlines()
            print("All output files:",ofnames)
            
            #find the log file and display log file
            cmd = "find %s -name *log.txt" % odir
            lines = os.popen(cmd).readlines()
            if len(lines) > 0:
                lfname = lines[0].rstrip()
                print("Log file is:", lfname)
                log = open(lfname).readlines()
                for line in log:         print(line.rstrip())
            
        else:
            print("Executable not found")
            
    
    In [2]:
    print("Marker : Calculation Started")
    compareSamples("data/bln69_minima_1.txt","data/bln69_energies_1.txt",\
                   "data/bln69_minima_2.txt","data/bln69_energies_2.txt")      
    print("Marker : Calculation Ended")
    
    Marker : Calculation Started
    Using executable /user/fcazals/home/projects/proj-soft/sbl-install/bin/sbl-emd-wp-euclid.exe
    
    All output files: ['sbl-wp-emd-euclid__emd_engine.xml\n', 'sbl-wp-emd-euclid__log.txt\n', 'sbl-wp-emd-euclid__transportation_plan.dot\n']
    Log file is: tmp-results-samples/sbl-wp-emd-euclid__log.txt
    Running:  sbl-emd-wp-euclid.exe --points-file data/bln69_minima_1.txt --points-file data/bln69_minima_2.txt -v -l -d tmp-results-samples --weights-file data/bln69_energies_1.txt --weights-file data/bln69_energies_2.txt
    
    Vertices Loader
    Statistics:
    Number of loaded data: 2
    -- Number of loaded weighted points in ensemble: 100
    -- Number of loaded weighted points in ensemble: 100
    
    EMD
    cumulative mass source: -9373.71124229750239464920
    cumulative mass demand: -9388.33313279270078055561
    STEP: dumping the mps file for the LP i.e. linear_program_solver_input.mps
    
    STEP: solving the LP with lp_solve...
    Done with exit code 512
    
    STEP: parsing the result from lp_solve i.e. linear_program_solver_input.res
    
    Statistics:
    EMD without connectivity constraint...
    
    Source flow: 0.00000 / -9373.71124 = -0.00000
    Remaining demand: 0.00000 / -9388.33313 = -0.00000
    Ratio flow/demand: -0.00000
    Transportation cost (normalized): 0.00000 / 0.00000 = -nan
    
    Report...
    
    End Run
    
    General Statistics:
    
    Times elapsed for computations (in seconds):
    -- EMD: 0.07592800000000000937
    Total: 0.07592800000000000937
    
    Marker : Calculation Ended
    

    Example : Comparing graphs

    Options.

    The main options of the compareGraphs method in the next cell are:

    • points1: Text file listing the first D-dimensional points
    • weights1: Text file listing the weights of the first points
    • edges1: Text file listing the edges between the first points
    • points2: Text file listing the second D-dimensional points
    • weights2: Text file listing the weights of the second points
    • edges2: Text file listing the edges between the second points
    • constraints: use EMD-CC
    • algorithm: with EMD-CC, select the algorithm (reg or star)
    • edgeSelection: with EMD-CC, selection scheme (min-cost or first-found)
    • recursion: with EMD-CC, recursion mode (norec, refined or coarse)
    In [3]:
    import re  #regular expressions
    import sys #misc system
    import os
    import pdb
    import shutil # python 3 only
    
    def compareGraphs(points1, weights1, edges1, points2, weights2, edges2,\
                      constraints = False, algorithm = "reg", edgeSelection = "min-cost",\
                      recursion = "refined"): 
        
        odir = "tmp-results-graph"
        if os.path.exists(odir):
            os.system("rm -rf %s" % odir)
        os.system( ("mkdir %s" % odir) )
        
        # check executable exists and is visible
        exe = shutil.which("sbl-emd-graph-euclid.exe")
        if exe:
            print(("Using executable %s\n" % exe))
            cmd = "sbl-emd-graph-euclid.exe --vertex-points-file %s --vertex-points-file %s -v -l -d %s" %\
                  (points1, points2, odir)
            if weights1 != None and weights2 != None:
                cmd += " --vertex-weights-file %s  --vertex-weights-file %s" % (weights1, weights2)
            if edges1 != None and edges2 != None:
                cmd += " --edges-file %s  --edges-file %s" % (edges1, edges2)
            if constraints:
                cmd += " --algorithm \"%s\" --edge-selection \"%s\" --recursion-mode \"%s\"" %\
                (algorithm, edgeSelection, recursion)
            os.system(cmd)
            
            cmd = "ls %s" % odir
            ofnames = os.popen(cmd).readlines()
            print("All output files:",ofnames)
            
            #find the log file and display log file
            cmd = "find %s -name *log.txt" % odir
            lines = os.popen(cmd).readlines()
            if len(lines) > 0:
                lfname = lines[0].rstrip()
                print("Log file is:", lfname)
                log = open(lfname).readlines()
                for line in log:         print(line.rstrip())
            
        else:
            print("Executable not found")
    
    In [4]:
    print("Marker : Calculation Started")
    compareGraphs("data/bln69_minima_1.txt","data/bln69_energies_1.txt", "data/bln69_edges_1.txt",\
                  "data/bln69_minima_2.txt","data/bln69_energies_2.txt", "data/bln69_edges_2.txt", True)      
    print("Marker : Calculation Ended")
    
    Marker : Calculation Started
    Using executable /user/fcazals/home/projects/proj-soft/sbl-install/bin/sbl-emd-graph-euclid.exe
    
    All output files: ['sbl-graph-emd-euclid__emd_engine.xml\n', 'sbl-graph-emd-euclid__log.txt\n', 'sbl-graph-emd-euclid__transportation_plan.dot\n']
    Log file is: tmp-results-graph/sbl-graph-emd-euclid__log.txt
    Running:  sbl-emd-graph-euclid.exe --vertex-points-file data/bln69_minima_1.txt --vertex-points-file data/bln69_minima_2.txt -v -l -d tmp-results-graph --vertex-weights-file data/bln69_energies_1.txt --vertex-weights-file data/bln69_energies_2.txt --edges-file data/bln69_edges_1.txt --edges-file data/bln69_edges_2.txt --algorithm reg --edge-selection min-cost --recursion-mode refined
    
    Vertices Loader
    Statistics:
    Number of loaded data: 2
    -- Number of loaded vertices / edges: 100 / 124
    -- Number of loaded vertices / edges: 100 / 123
    
    EMD
    cumulative mass source: -9373.71124229750239464920
    cumulative mass demand: -9388.33313279270078055561
    STEP: dumping the mps file for the LP i.e. linear_program_solver_input.mps
    
    STEP: solving the LP with lp_solve...
    Done with exit code 512
    
    STEP: parsing the result from lp_solve i.e. linear_program_solver_input.res
    
    Statistics:
    EMD without connectivity constraint...
    
    Source flow: 0.00000 / -9373.71124 = -0.00000
    Remaining demand: 0.00000 / -9388.33313 = -0.00000
    Ratio flow/demand: -0.00000
    Transportation cost (normalized): 0.00000 / 0.00000 = -nan
    
    Report...
    
    End Run
    
    General Statistics:
    
    Times elapsed for computations (in seconds):
    -- EMD: 0.08190300000000000358
    Total: 0.08190300000000000358
    
    Marker : Calculation Ended