 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 >> )  . Such problems were later rephrased in terms of Riemannian geometry and measure theory  , 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  , statistics  , geometry processing  , biophysics  ...

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 , and the demand nodes respectively. We assume that these nodes are embedded in some metric space , and denote a metric in this space. We assume that each source node has a mass of supply (aka weights) denoted , and equivalently, that each demand node has a mass of demand . The sums of masses for all source and demand nodes are denoted and , that is and . We also assume that transporting a unit mass of supply from to demand node has unit cost .

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 to is defined by triples , with the so-called flow quantities. Note that there are at most such triples.

Finding the optimal i.e. least cost transport plan amounts to solving the following linear program (LP): 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  : ## 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 and . 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 to export flow to a connected component of . 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 from . Let be the subgraph of induced by the vertices from such that for each such vertex , there is at least one edge emanating from a vertex of the subgraph with . The transport plan is called valid iff the subgraph 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 induces, through edges carrying strictly positive flow, a connected set of vertices from .

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  . This motivates the definition of a new problem denoted .

The problem solved: .

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 ( ). 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 . We define the following upper bound for the maximum total cost of any admissible flow as Consider an upper bound , and an upper bound on the number of edges which carry non null flow. Denote and the set of all connected subgraphs of and , respectively. Problem is defined as follows: 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  . Practically, various solvers can be used, e.g. the one from the Computational Geometry Algorithms Library  , 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 .

## Earth Mover Distance with connectivity constraints

Algorithms.

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

• the decision version of is NP-complete,

• a Polynomial Time Approximation Scheme for exists when transport plans involving a large number of edges are allowed.

Practically, we solve with two algorithms denoted and . 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 . A greedy strategy returning an admissible solution, given a maximum cost provided by the user.

• Recursion mode refined: algorithm . Given a range 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 . 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: It can be shown that the Earth Mover Distance as defined by Eq. (emd-eq-emd-lp) yields a metric, provided that is itself a metric, and that the sum of weights for the source and the demand graphs are equal  . On the other hand, earth mover distances with connectivity constraints are not symmetric in general  , 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: # 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 :
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.rstrip()
print("Log file is:", lfname)
log = open(lfname).readlines()
for line in log:         print(line.rstrip())

else:
print("Executable not found")


In :
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 :
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.rstrip()
print("Log file is:", lfname)
log = open(lfname).readlines()
for line in log:         print(line.rstrip())

else:
print("Executable not found")

In :
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