From 2ed9afc052bee7956f6abb195947de1f80cb9d91 Mon Sep 17 00:00:00 2001
From: Arnur Nigmetov <anigmetov@lbl.gov>
Date: Fri, 7 May 2021 22:38:27 -0700
Subject: Fix bugs in matching distance for modules.

1. Twist reduction should not be used for modules, clearing makes no
   sense and creates erroneous points. Just use standard reduction.
2. Infinite points for modules (gens/relations case) are treated
   slightly differently.
---
 matching/include/matching_distance.hpp            |  1 +
 matching/include/persistence_module.h             |  1 +
 matching/include/persistence_module.hpp           |  2 +-
 matching/include/phat/compute_persistence_pairs.h | 19 ++++++++++++++-----
 4 files changed, 17 insertions(+), 6 deletions(-)

diff --git a/matching/include/matching_distance.hpp b/matching/include/matching_distance.hpp
index f7f44a5..9beab1f 100644
--- a/matching/include/matching_distance.hpp
+++ b/matching/include/matching_distance.hpp
@@ -362,6 +362,7 @@ namespace md {
         // TODO: think about this - how to call Hera
         auto dgm_a = module_a_.weighted_slice_diagram(line);
         auto dgm_b = module_b_.weighted_slice_diagram(line);
+
         R result;
         if (params_.hera_epsilon > static_cast<R>(0)) {
             result = hera::bottleneckDistApprox(dgm_a, dgm_b, params_.hera_epsilon) / ( params_.hera_epsilon + 1);
diff --git a/matching/include/persistence_module.h b/matching/include/persistence_module.h
index b68c21e..4df2148 100644
--- a/matching/include/persistence_module.h
+++ b/matching/include/persistence_module.h
@@ -11,6 +11,7 @@
 
 #include "phat/boundary_matrix.h"
 #include "phat/compute_persistence_pairs.h"
+#include "phat/algorithms/standard_reduction.h"
 
 #include "common_util.h"
 #include "dual_point.h"
diff --git a/matching/include/persistence_module.hpp b/matching/include/persistence_module.hpp
index 7479e02..022b43d 100644
--- a/matching/include/persistence_module.hpp
+++ b/matching/include/persistence_module.hpp
@@ -158,7 +158,7 @@ namespace md {
         get_slice_projection_matrix(slice, phat_matrix, gen_projections, rel_projections);
 
         phat::persistence_pairs phat_persistence_pairs;
-        phat::compute_persistence_pairs<phat::twist_reduction>(phat_persistence_pairs, phat_matrix);
+        phat::compute_persistence_pairs<phat::standard_reduction>(phat_persistence_pairs, phat_matrix, true);
 
         Diagram<Real> dgm;
 
diff --git a/matching/include/phat/compute_persistence_pairs.h b/matching/include/phat/compute_persistence_pairs.h
index 06f5372..bd5e3ac 100644
--- a/matching/include/phat/compute_persistence_pairs.h
+++ b/matching/include/phat/compute_persistence_pairs.h
@@ -93,7 +93,7 @@ namespace phat {
     }
 
     template< typename ReductionAlgorithm, typename Representation >
-    void compute_persistence_pairs( persistence_pairs& pairs, boundary_matrix< Representation >& boundary_matrix ) {
+    void compute_persistence_pairs( persistence_pairs& pairs, boundary_matrix< Representation >& boundary_matrix, bool for_persistence_module = false ) {
         ReductionAlgorithm reduce;
         reduce( boundary_matrix );
         pairs.clear();
@@ -107,10 +107,19 @@ namespace phat {
                 pairs.append_pair( birth, death );
             }
         }
-        // infinite pairs: column idx is 0, and row idx does not contain a lowest one
-        for( index idx = 0; idx < boundary_matrix.get_num_cols(); idx++ ) {
-            if(boundary_matrix.is_empty(idx) && max_indices.count(idx) == 0 ) {
-                pairs.append_pair( idx, k_infinity_index);
+        if (for_persistence_module) {
+             // infinite pairs for modules: if generator is not matched to relation, it lives forever
+            for( index idx = 0; idx < boundary_matrix.get_num_rows(0); idx++ ) {
+                if(max_indices.count(idx) == 0 ) {
+                    pairs.append_pair( idx, k_infinity_index);
+                }
+            }
+        } else {
+            // infinite pairs: column idx is 0, and row idx does not contain a lowest one
+            for( index idx = 0; idx < boundary_matrix.get_num_cols(); idx++ ) {
+                if(boundary_matrix.is_empty(idx) && max_indices.count(idx) == 0 ) {
+                    pairs.append_pair( idx, k_infinity_index);
+                }
             }
         }
      }
-- 
cgit v1.2.3