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 >> ) [130] . Such problems were later rephrased in terms of Riemannian geometry and measure theory [170] , 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 [154] , statistics [117] , geometry processing [163] , biophysics [45] ...
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[154] :
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 [39] . 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 [105] . Practically, various solvers can be used, e.g. the one from the Computational Geometry Algorithms Library [50] , lp_solve, the CPLEX solver from IBM, etc.
This package currently proposes two options to solve the LP problems of Eq. emd-eq-emd-lp :
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 [154] . On the other hand, earth mover distances with connectivity constraints are not symmetric in general [39] , 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:
SBL::CADS::T_Earth_mover_distance_vertices_accessor_vector < VertexType , MassType > for representing the source or the demand by a STL vector of objects. The object type is defined by the template argument VertexType while MassType defines the mass associated to each vertex.
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.
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.
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.
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.
Earth Mover Distance module
The following example shows a simple instantiation of the EMD module within a workflow.
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:
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).
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
Main option for the LP program used by EMD. One needs to specify the linear solver used:
–lp-solverstring: Linear solver used; options are: lp_solve or clp
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-filestring: Description of points in point d format – one point per line –weights-filestring: Description of positive weights – the ordering as points –edges-filestring: File listing edges of the graph – two integers in 0..n-1
Main options for EMD with connectivity constraints. The main options are:
–constraintsstring: use EMD-CC –algorithmstring: with EMD-CC –select the algorithm (reg or star) –edgeSelectionstring: with EMD-CC – selection scheme (min-cost or first-found) –recursionstring: with EMD-CC – recursion mode (norec / refined / coarse)
Sum of values:0.9999999999999998
Sum of values:0.9999999999999998
Using executable /user/fcazals/home/projects/proj-soft/sbl-install/bin/sbl-emd-wp-lrmsd.exe
Executing /user/fcazals/home/projects/proj-soft/sbl-install/bin/sbl-emd-wp-lrmsd.exe --points-file ../demos/data/bln69_minima_1.txt --points-file ../demos/data/bln69_minima_2.txt -v -l -d results-samples-bln69-lp_solve --weights-file bln69_energies_1-weights.txt --weights-file bln69_energies_2-weights.txt --lp-solver lp_solve
All output files: ['sbl-emd-wp-lrmsd__emd_engine.xml\n', 'sbl-emd-wp-lrmsd__log.txt\n', 'sbl-emd-wp-lrmsd__transportation_plan.dot\n']
Showing file results-samples-bln69-lp_solve/sbl-emd-wp-lrmsd__log.txt
++Showing file results-samples-bln69-lp_solve/sbl-emd-wp-lrmsd__log.txt
Running: /user/fcazals/home/projects/proj-soft/sbl-install/bin/sbl-emd-wp-lrmsd.exe --points-file ../demos/data/bln69_minima_1.txt --points-file ../demos/data/bln69_minima_2.txt -v -l -d results-samples-bln69-lp_solve --weights-file bln69_energies_1-weights.txt --weights-file bln69_energies_2-weights.txt --lp-solver lp_solve
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
Starting module SBL::Modules::T_Earth_mover_distance_module...
cumulative mass source: 0.99999999999999977796
cumulative mass demand: 0.99999999999999977796
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 0
STEP: parsing the result from lp_solve i.e. linear_program_solver_input.res
Statistics:
EMD without connectivity constraint...
Source flow: 1.00000 / 1.00000 = 1.00000
Remaining demand: 0.00000 / 1.00000 = 0.00000
Ratio flow/demand: 1.00000
Transportation cost (normalized): 0.16320 / 1.00000 = 0.16320
Report...
Ending module SBL::Modules::T_Earth_mover_distance_module...
End Run
General Statistics:
Times elapsed for computations (in seconds):
-- EMD: 1.11253099999999993663
Total: 1.11253099999999993663
--Done
XML: 1 / 1 files were loaded
EMD_comparators: summary
Total cost: 0.16319675899339775
Total flow: 0.9999998508942526
Normalized cost: 0.1631967833269761
Using executable /user/fcazals/home/projects/proj-soft/sbl-install/bin/sbl-emd-wp-lrmsd.exe
Executing /user/fcazals/home/projects/proj-soft/sbl-install/bin/sbl-emd-wp-lrmsd.exe --points-file ../demos/data/bln69_minima_1.txt --points-file ../demos/data/bln69_minima_2.txt -v -l -d results-samples-bln69-clp --weights-file bln69_energies_1-weights.txt --weights-file bln69_energies_2-weights.txt --lp-solver clp
All output files: ['sbl-emd-wp-lrmsd__emd_engine.xml\n', 'sbl-emd-wp-lrmsd__log.txt\n', 'sbl-emd-wp-lrmsd__transportation_plan.dot\n']
Showing file results-samples-bln69-clp/sbl-emd-wp-lrmsd__log.txt
++Showing file results-samples-bln69-clp/sbl-emd-wp-lrmsd__log.txt
Running: /user/fcazals/home/projects/proj-soft/sbl-install/bin/sbl-emd-wp-lrmsd.exe --points-file ../demos/data/bln69_minima_1.txt --points-file ../demos/data/bln69_minima_2.txt -v -l -d results-samples-bln69-clp --weights-file bln69_energies_1-weights.txt --weights-file bln69_energies_2-weights.txt --lp-solver clp
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
Starting module SBL::Modules::T_Earth_mover_distance_module...
cumulative mass source: 0.99999999999999977796
cumulative mass demand: 0.99999999999999977796
STEP: dumping the mps file for the LP i.e. linear_program_solver_input.mps
STEP: solving the LP with clp...
Done with exit code 0
STEP: parsing the result from clp i.e. linear_program_solver_input.txt
Statistics:
EMD without connectivity constraint...
Source flow: 0.99886 / 1.00000 = 0.99886
Remaining demand: 0.00114 / 1.00000 = 0.00114
Ratio flow/demand: 0.99886
Transportation cost (normalized): 0.16262 / 0.99886 = 0.16280
Report...
Ending module SBL::Modules::T_Earth_mover_distance_module...
End Run
General Statistics:
Times elapsed for computations (in seconds):
-- EMD: 1.10677999999999987502
Total: 1.10677999999999987502
--Done
XML: 1 / 1 files were loaded
EMD_comparators: summary
Total cost: 0.16261716647313668
Total flow: 0.9988600136857713
Normalized cost: 0.1628027593907608