Fix multiple bugs in matching distance for modules, add example.

7 files changed, 89 insertions, 14 deletions

diff --git a/matching/CMakeLists.txt b/matching/CMakeLists.txt
index 3ee0f6b..a391d84 100644
--- a/matching/CMakeLists.txt
+++ b/matching/CMakeLists.txt
@@ -48,5 +48,9 @@ endif()
 add_executable(matching_dist "example/matching_dist.cpp" ${MD_HEADERS} ${BT_HEADERS} )
 target_link_libraries(matching_dist PUBLIC ${libraries})
 
+add_executable(module_example "example/module_example.cpp" ${MD_HEADERS} ${BT_HEADERS} )
+target_link_libraries(module_example PUBLIC ${libraries})
+
+
 add_executable(matching_distance_test ${SRC_TEST_FILES} ${BT_HEADERS} ${MD_HEADERS})
 target_link_libraries(matching_distance_test PUBLIC ${libraries})
diff --git a/matching/README.md b/matching/README.md
index 7dc1874..fc97441 100644
--- a/matching/README.md
+++ b/matching/README.md
@@ -1,5 +1,69 @@
-Matching distance between bifiltrations.
+# Matching distance between bifiltrations and 2-persistence modules.
 
-Currently supports only 1-critical bi-filtrations
-in PHAT-like format: boundary matrix + critical values 
-of a simplex in each row
+## Accompanying paper
+M. Kerber, A. Nigmetov, 
+Efficient Approximation of the Matching Distance for 2-parameter persistence.
+SoCG 2020.
+
+Bug reports can be sent to "anigmetov EMAIL SIGN lbl DOT gov".
+
+## Dependencies
+
+* Your compiler must support C++11.
+* Boost.
+
+## Usage:
+
+1. To use a standalone command-line utility matching_dist:
+
+`./matching_dist -d dimension -e relative_error file1 file2`
+
+If `relative_error` is not specified, the default 0.1 is used.
+If `dimension` is not specified, the default value is 0.
+Run `./matching_dist` without parameters to see other options.
+
+The output is an approximation of the exact distance (which is assumed to be non-zero).
+Precisely: if  *d_exact* is the true distance and *d_approx* is the output, then 
+
+>     | d_exact - d_approx | / d_exact < relative_error.
+
+Files file1 and file2 must contain 1-critical bi-filtrations in a plain text format which is similar to PHAT. The first line of a file must say *bifiltration_phat_like*. The second line contains the total number of simplices *N*. The next *N* lines contain simplices in the format *dim x y boundary*.
+* *dim*: the dimension of the simplex
+* *x, y*: coordinates of the critical value
+* *boundary*: indices of simplices forming the boundary of the current simplex. Indices are separated by space.
+* Simplices are indexed starting from 0.
+
+For example, the bi-filtration of a segment with vertices appearing at (0,0) and the 1-segment appearing at (3,4) shall be written as:
+
+>    bifiltration_phat_like
+>    3
+>    \# lines starting with \# are ignored
+>    \# vertex A has dimension 0, hence no boundary, its index is 0
+>    0 0 0
+>    \# vertex B has index 1
+>    0 0 0
+>    \# 1-dimensional simplex {A, B}
+>    1 3 4 0 1
+
+2. To use from your code.
+
+Here you can compute the matching distance either between bi-filtrations or between  persistence modules.
+First, you need to include `#include "matching_distance.h"` Practically every class you need is parameterized by Real type, which should be either float or double.  The header provides two functions called `matching_distance.` 
+See `example/module_example.cpp` for additional details.
+
+## License
+
+See `licence.txt` in the repository root folder.
+
+## Building
+
+CMakeLists.txt can be used to build the command-line utility in the standard
+way.  On Linux/Mac/Windows with Cygwin:
+>    `mkdir build`
+>    `cd build`
+>   `cmake ..`
+>    `make`
+
+On Windows with Visual Studio: use `cmake-gui` to create the solution in build directory and build it with VS.
+
+The library itself is header-only and does not require separate compilation.
diff --git a/matching/example/matching_dist.cpp b/matching/example/matching_dist.cpp
index 13e5c6d..734f6cc 100644
--- a/matching/example/matching_dist.cpp
+++ b/matching/example/matching_dist.cpp
@@ -13,9 +13,6 @@
 #include "spdlog/spdlog.h"

 #include "spdlog/fmt/ostr.h"

 

-//#include "persistence_module.h"

-#include "bifiltration.h"

-#include "box.h"

 #include "matching_distance.h"

 

 using Real = double;

diff --git a/matching/include/matching_distance.h b/matching/include/matching_distance.h
index 276cc3a..e82a97c 100644
--- a/matching/include/matching_distance.h
+++ b/matching/include/matching_distance.h
@@ -152,7 +152,7 @@ namespace md {
         int max_depth {6}; // maximal number of refinenemnts
         int initialization_depth {3};
         int dim {0}; // in which dim to calculate the distance; use ALL_DIMENSIONS to get max over all dims
-        BoundStrategy bound_strategy {BoundStrategy::bruteforce};
+        BoundStrategy bound_strategy {BoundStrategy::local_combined};
         TraverseStrategy traverse_strategy {TraverseStrategy::breadth_first};
         bool tolerate_max_iter_exceeded {true};
         Real actual_error {std::numeric_limits<Real>::max()};
diff --git a/matching/include/matching_distance.hpp b/matching/include/matching_distance.hpp
index 48c8464..7cff073 100644
--- a/matching/include/matching_distance.hpp
+++ b/matching/include/matching_distance.hpp
@@ -193,6 +193,7 @@ namespace md {
         // make all coordinates non-negative
         auto min_coord = std::min(module_a_.minimal_coordinate(),
                 module_b_.minimal_coordinate());
+        spd::debug("in DistanceCalculator ctor, min_coord = {}", min_coord);
         if (min_coord < 0) {
             module_a_.translate(-min_coord);
             module_b_.translate(-min_coord);
diff --git a/matching/include/persistence_module.h b/matching/include/persistence_module.h
index e99771f..4a261bb 100644
--- a/matching/include/persistence_module.h
+++ b/matching/include/persistence_module.h
@@ -47,7 +47,7 @@ namespace md {
             IndexVec components_;
 
             Relation() {}
-            Relation(const Point<Real>& _pos, const IndexVec& _components);
+            Relation(const Point<Real>& _pos, const IndexVec& _components) : position_(_pos), components_(_components) {}
 
             Real get_x() const { return position_.x; }
             Real get_y() const { return position_.y; }
diff --git a/matching/include/persistence_module.hpp b/matching/include/persistence_module.hpp
index 6e49b2e..128fed9 100644
--- a/matching/include/persistence_module.hpp
+++ b/matching/include/persistence_module.hpp
@@ -34,10 +34,10 @@ namespace md {
     template<class Real>
     void ModulePresentation<Real>::init_boundaries()
     {
-        max_x_ = std::numeric_limits<Real>::max();
-        max_y_ = std::numeric_limits<Real>::max();
-        min_x_ = -std::numeric_limits<Real>::max();
-        min_y_ = -std::numeric_limits<Real>::max();
+        max_x_ = -std::numeric_limits<Real>::max();
+        max_y_ = -std::numeric_limits<Real>::max();
+        min_x_ = std::numeric_limits<Real>::max();
+        min_y_ = std::numeric_limits<Real>::max();
 
         for(const auto& gen : positions_) {
             min_x_ = std::min(gen.x, min_x_);
@@ -55,6 +55,7 @@ namespace md {
         generators_(_generators),
         relations_(_relations)
     {
+        positions_ = concat_gen_and_rel_positions(generators_, relations_);
         init_boundaries();
     }
 
@@ -85,6 +86,7 @@ namespace md {
     void ModulePresentation<Real>::project_generators(const DualPoint<Real>& slice,
             IndexVec& sorted_indices, RealVec& projections) const
     {
+        spd::debug("Enter project_generators, slice = {}", slice);
         size_t num_gens = generators_.size();
 
         RealVec gen_values;
@@ -97,6 +99,7 @@ namespace md {
         projections.reserve(num_gens);
         for(auto i : sorted_indices) {
             projections.push_back(gen_values[i]);
+            spd::debug("added push = {}", gen_values[i]);
         }
     }
 
@@ -104,6 +107,8 @@ namespace md {
     void ModulePresentation<Real>::project_relations(const DualPoint<Real>& slice, IndexVec& sorted_rel_indices,
             RealVec& projections) const
     {
+
+        spd::debug("Enter project_relations, slice = {}", slice);
         size_t num_rels = relations_.size();
 
         RealVec rel_values;
@@ -116,12 +121,14 @@ namespace md {
         projections.reserve(num_rels);
         for(auto i : sorted_rel_indices) {
             projections.push_back(rel_values[i]);
+            spd::debug("added push = {}", rel_values[i]);
         }
     }
 
     template<class Real>
     Diagram<Real> ModulePresentation<Real>::weighted_slice_diagram(const DualPoint<Real>& slice) const
     {
+        spd::debug("Enter weighted_slice_diagram, slice = {}", slice);
         IndexVec sorted_gen_indices, sorted_rel_indices;
         RealVec gen_projections, rel_projections;
 
@@ -138,7 +145,8 @@ namespace md {
                 j = sorted_gen_indices[j];
             }
             std::sort(current_relation.begin(), current_relation.end());
-            phat_matrix.set_dim(i, current_relation.size());
+            // modules do not have dimension, set all to 0
+            phat_matrix.set_dim(i, 0);
             phat_matrix.set_col(i, current_relation);
         }
 
@@ -154,6 +162,7 @@ namespace md {
             bool is_finite_pair = new_pair.second != phat::k_infinity_index;
             Real birth = gen_projections.at(new_pair.first);
             Real death = is_finite_pair ? rel_projections.at(new_pair.second) : real_inf;
+            spd::debug("i = {}, birth = {}, death = {}", i, new_pair.first, new_pair.second);
             if (birth != death) {
                 dgm.emplace_back(birth, death);
             }