store diameters in reduction matrix

2 files changed, 76 insertions, 66 deletions

diff --git a/ripser.cpp b/ripser.cpp
index 902caed..7fbbe66 100644
--- a/ripser.cpp
+++ b/ripser.cpp
@@ -193,43 +193,62 @@ inline void set_coefficient(index_t& e, const coefficient_t c) { e = c; }
 
 inline const entry_t& get_entry(const entry_t& e) { return e; }
 
+template <typename Entry>
 struct greater_index {
-    bool operator() (const entry_t& a, const entry_t& b) {
+    bool operator() (const Entry& a, const Entry& b) {
         return get_index(a) > get_index(b);
     }
 };
 
 
 #ifdef STORE_DIAMETERS
-class entry_diameter_t: public std::pair<entry_t, value_t> {
+typedef std::pair<value_t, index_t> diameter_index_t;
+inline value_t get_diameter(diameter_index_t i) {
+    return i.first;
+}
+inline index_t get_index(diameter_index_t i) {
+    return i.second;
+}
+#else
+typedef index_t diameter_index_t;
+#endif
+
+#ifdef STORE_DIAMETERS
+class diameter_entry_t: public std::pair<value_t, entry_t> {
 public:
-    entry_diameter_t(std::pair<entry_t, value_t> p) : std::pair<entry_t, value_t>(p) {}
+    diameter_entry_t(std::pair<value_t, entry_t> p) : std::pair<value_t, entry_t>(p) {}
 #ifdef USE_COEFFICIENTS
-    entry_diameter_t(entry_t e) : std::pair<entry_t, value_t>(e, 0) {}
+    diameter_entry_t(entry_t e) : std::pair<value_t, entry_t>(0, e) {}
 #endif
-    entry_diameter_t(index_t i) : std::pair<entry_t, value_t>(i, 0) {}
+    diameter_entry_t(index_t i) : std::pair<value_t, entry_t>(0, i) {}
 };
 
-inline const entry_t& get_entry(const entry_diameter_t& p) { return p.first; }
-inline entry_t& get_entry(entry_diameter_t& p) { return p.first; }
-inline const index_t get_index(const entry_diameter_t& p) { return get_index(get_entry(p)); }
-inline const coefficient_t get_coefficient(const entry_diameter_t& p) { return get_coefficient(get_entry(p)); }
-inline const value_t& get_diameter(const entry_diameter_t& p) { return p.second; }
-inline void set_coefficient(entry_diameter_t& p, const coefficient_t c) { set_coefficient(get_entry(p), c); }
-inline entry_diameter_t make_entry_diameter(index_t _index, coefficient_t _coefficient, value_t _diameter) {
-    return std::make_pair(make_entry(_index, _coefficient), _diameter);
+inline const entry_t& get_entry(const diameter_entry_t& p) { return p.second; }
+inline entry_t& get_entry(diameter_entry_t& p) { return p.second; }
+inline const index_t get_index(const diameter_entry_t& p) { return get_index(get_entry(p)); }
+inline const coefficient_t get_coefficient(const diameter_entry_t& p) { return get_coefficient(get_entry(p)); }
+inline const value_t& get_diameter(const diameter_entry_t& p) { return p.first; }
+inline void set_coefficient(diameter_entry_t& p, const coefficient_t c) { set_coefficient(get_entry(p), c); }
+inline diameter_entry_t make_diameter_entry(value_t _diameter, index_t _index, coefficient_t _coefficient) {
+    return std::make_pair(_diameter, make_entry(_index, _coefficient));
+}
+inline diameter_entry_t make_diameter_entry(diameter_index_t _diameter_index, coefficient_t _coefficient) {
+    return std::make_pair(get_diameter(_diameter_index), make_entry(get_index(_diameter_index), _coefficient));
 }
 
 
 struct greater_diameter_or_index {
-    inline bool operator() (const entry_diameter_t& a, const entry_diameter_t& b) {
+    inline bool operator() (const diameter_entry_t& a, const diameter_entry_t& b) {
         return (get_diameter(a) > get_diameter(b)) ||
             ((get_diameter(a) == get_diameter(b)) && (get_index(a) > get_index(b)));
     }
 };
 #else
-typedef entry_t entry_diameter_t;
-inline entry_diameter_t make_entry_diameter(index_t _index, coefficient_t _coefficient, value_t _diameter) {
+typedef entry_t diameter_entry_t;
+inline diameter_entry_t make_diameter_entry(value_t _diameter, index_t _index, coefficient_t _coefficient) {
+    return make_entry(_index, _coefficient);
+}
+inline diameter_entry_t make_diameter_entry(index_t _index, coefficient_t _coefficient) {
     return make_entry(_index, _coefficient);
 }
 #endif
@@ -333,8 +352,9 @@ public:
         }
         auto result = std::make_pair(make_entry(
             modified_idx + binomial_coeff( n , k + 1 ),
-            k & 1 ? -1 : 1
-        ), n);
+            k & 1 ? -1 : 1),
+            n);
+        
         --n;
         return result;
     }
@@ -552,7 +572,7 @@ public:
 
 #ifdef USE_COEFFICIENTS
 template <typename Heap>
-inline entry_diameter_t pop_pivot(Heap& column, coefficient_t modulus)
+inline diameter_entry_t pop_pivot(Heap& column, coefficient_t modulus)
 {
     
     if( column.empty() )
@@ -583,7 +603,7 @@ inline entry_diameter_t pop_pivot(Heap& column, coefficient_t modulus)
 #else
 
 template <typename Heap>
-inline entry_diameter_t pop_pivot(Heap& column, coefficient_t modulus)
+inline diameter_entry_t pop_pivot(Heap& column, coefficient_t modulus)
 {
     
     if( column.empty() )
@@ -607,9 +627,9 @@ inline entry_diameter_t pop_pivot(Heap& column, coefficient_t modulus)
 #endif
 
 template <typename Heap>
-inline entry_diameter_t get_pivot(Heap& column, coefficient_t modulus)
+inline diameter_entry_t get_pivot(Heap& column, coefficient_t modulus)
 {
-    entry_diameter_t result = pop_pivot(column, modulus);
+    diameter_entry_t result = pop_pivot(column, modulus);
     if( get_index(result) != -1 ) {
         column.push(result);
     }
@@ -617,17 +637,6 @@ inline entry_diameter_t get_pivot(Heap& column, coefficient_t modulus)
 }
 
 
-#ifdef STORE_DIAMETERS
-typedef std::pair<value_t, index_t> diameter_index_t;
-inline value_t get_diameter(diameter_index_t i) {
-    return i.first;
-}
-inline index_t get_index(diameter_index_t i) {
-    return i.second;
-}
-#else
-typedef index_t diameter_index_t;
-#endif
 
 template <typename Comparator>
 void assemble_columns_to_reduce (
@@ -745,7 +754,7 @@ template <typename Heap>
 inline void push_entry(Heap& column, index_t i, coefficient_t c, value_t diameter) {
     entry_t e = make_entry(i, c);
 #ifdef STORE_DIAMETERS
-    column.push(std::make_pair(e, diameter));
+    column.push(std::make_pair(diameter, e));
 #else
     column.push(e);
 #endif
@@ -769,27 +778,27 @@ void compute_pairs(
     #endif
     
     #ifdef ASSEMBLE_REDUCTION_MATRIX
-    compressed_sparse_matrix <entry_t> reduction_matrix;
+    compressed_sparse_matrix <diameter_entry_t> reduction_matrix;
     #else
     #ifdef USE_COEFFICIENTS
-    std::vector <entry_t> reduction_coefficients;
+    std::vector <diameter_entry_t> reduction_coefficients;
     #endif
     #endif
     
 //    size_t boundary_additions = 0;
 
     for (index_t i = 0; i < columns_to_reduce.size(); ++i) {
-        index_t column_to_reduce = get_index(columns_to_reduce[i]);
+        auto column_to_reduce = columns_to_reduce[i];
         
         #ifdef ASSEMBLE_REDUCTION_MATRIX
         
-        std::priority_queue<entry_t, std::vector<entry_t>, greater_index> reduction_column;
+        std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>, greater_index<diameter_entry_t>> reduction_column;
         
         #endif
         
         
         #ifdef STORE_DIAMETERS
-        std::priority_queue<entry_diameter_t, std::vector<entry_diameter_t>, greater_diameter_or_index > working_coboundary;
+        std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>, greater_diameter_or_index > working_coboundary;
         #else
         std::priority_queue<entry_t, std::vector<entry_t>, decltype(comp) > working_coboundary(comp);
         #endif
@@ -797,7 +806,7 @@ void compute_pairs(
         #ifdef STORE_DIAMETERS
         value_t diameter = get_diameter(columns_to_reduce[i]);
         #else
-        value_t diameter = comp_prev.diameter(column_to_reduce);
+        value_t diameter = comp_prev.diameter(get_index(column_to_reduce));
         #endif
 
         
@@ -809,13 +818,13 @@ void compute_pairs(
         
         index_t j = i;
         
-        index_t column_to_add = column_to_reduce;
+        auto column_to_add = column_to_reduce;
         
         
         // start with a dummy pivot entry with coefficient -1 in order to initialize working_coboundary
         // with the coboundary of the simplex with index column_to_reduce
         
-        entry_diameter_t pivot = make_entry_diameter(-1, modulus - 1, 0);
+        diameter_entry_t pivot = make_diameter_entry(0, -1, modulus - 1);
         
 //        std::cout << "reducing " << column_to_reduce << ":" << std::endl;
         
@@ -826,10 +835,10 @@ void compute_pairs(
 
         // initialize reduction_matrix as identity matrix
         #ifdef ASSEMBLE_REDUCTION_MATRIX
-        reduction_matrix.push_back(make_entry(column_to_reduce, 1));
+        reduction_matrix.push_back(make_diameter_entry(column_to_reduce, 1));
         #else
         #ifdef USE_COEFFICIENTS
-        reduction_coefficients.push_back(make_entry(column_to_reduce, 1));
+        reduction_coefficients.push_back(make_diameter_entry(column_to_reduce, 1));
         #endif
         #endif
         
@@ -840,7 +849,7 @@ void compute_pairs(
             const coefficient_t factor = modulus - get_coefficient(pivot);
 
 
-//            std::priority_queue<entry_diameter_t, std::vector<entry_diameter_t>, decltype(comp) > eliminating_coboundary(comp);
+//            std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>, decltype(comp) > eliminating_coboundary(comp);
 
 //            std::cout << "w:" << get_column_vector(working_coboundary, modulus) << std::endl;
 
@@ -853,22 +862,28 @@ void compute_pairs(
 
 
                 #ifdef ASSEMBLE_REDUCTION_MATRIX
-                entry_t simplex = *it;
+                auto simplex = *it;
                 coefficient_t simplex_coefficient = get_coefficient(simplex);
                 simplex_coefficient *= factor;
                 simplex_coefficient %= modulus;
-                reduction_column.push( make_entry( get_index(simplex), simplex_coefficient ) );
                 #else
                 #ifdef USE_COEFFICIENTS
-                const entry_t& simplex = reduction_coefficients[j];
+                const auto& simplex = reduction_coefficients[j];
                 #else
-                const entry_t simplex = column_to_add;
+                const auto simplex = column_to_add;
                 #endif
                 #endif
                 
-                //TODO: cache instead of compute?
+                #ifdef STORE_DIAMETERS
+                value_t simplex_diameter = get_diameter(simplex);
+                #else
                 value_t simplex_diameter = comp_prev.diameter(get_index(simplex));
+                #endif
+                assert(simplex_diameter == comp_prev.diameter(get_index(simplex)));
 
+                #ifdef ASSEMBLE_REDUCTION_MATRIX
+                reduction_column.push( make_diameter_entry( simplex_diameter, get_index(simplex), simplex_coefficient ) );
+                #endif
                 
                 std::vector<index_t> vertices = vertices_of_simplex(get_index(simplex), dim, n, binomial_coeff);
                 
@@ -878,27 +893,20 @@ void compute_pairs(
                     
                     entry_t coface = coface_descriptor.first;
                     index_t covertex = coface_descriptor.second;
-//
+
                     index_t coface_index = get_index(coface);
-//
+
                     value_t coface_diameter = simplex_diameter;
                     for (index_t v : vertices) {
                         coface_diameter = std::max(coface_diameter, dist(v, covertex));
                     }
                     
-//                    value_t coface_diameter = comp.diameter(coface_index);
-                    
-//                    assert(coface_diameter_variant == coface_diameter);
+                    assert(comp.diameter(coface_index) == coface_diameter);
                     
                     if (coface_diameter <= threshold) {
-                        coefficient_t coface_coefficient = get_coefficient(coface) + modulus;
-                        coface_coefficient %= modulus;
-                        
-                        coface_coefficient *= get_coefficient(simplex);
-                        coface_coefficient %= modulus;
-                        
-                        coface_coefficient *= factor;
+                        coefficient_t coface_coefficient = get_coefficient(coface) * get_coefficient(simplex) * factor;
                         coface_coefficient %= modulus;
+                        if (coface_coefficient < 0) coface_coefficient += modulus;
                         
                         assert(coface_coefficient >= 0);
                         
@@ -941,18 +949,18 @@ void compute_pairs(
                     #ifdef ASSEMBLE_REDUCTION_MATRIX
                     reduction_matrix.pop_back();
                     while (true) {
-                        entry_t e = get_entry(pop_pivot(reduction_column, modulus));
+                        diameter_entry_t e = pop_pivot(reduction_column, modulus);
                         index_t index = get_index(e);
                         if (index == -1)
                             break;
                         const coefficient_t coefficient = inverse * get_coefficient(e) % modulus;
                         assert(coefficient > 0);
-                        reduction_matrix.push_back(make_entry(index, coefficient));
+                        reduction_matrix.push_back(make_diameter_entry(get_diameter(e), index, coefficient));
                     }
                     #else
                     #ifdef USE_COEFFICIENTS
                     reduction_coefficients.pop_back();
-                    reduction_coefficients.push_back(make_entry(column_to_reduce, inverse));
+                    reduction_coefficients.push_back(make_diameter_entry(column_to_reduce, inverse));
                     #endif
                     #endif
                     
@@ -977,7 +985,7 @@ void compute_pairs(
                 }
                 
                 j = pair->second;
-                column_to_add = get_index(columns_to_reduce[j]);
+                column_to_add = columns_to_reduce[j];
             }
 
         } while ( get_index(pivot) != -1 );
diff --git a/ripser.xcodeproj/project.pbxproj b/ripser.xcodeproj/project.pbxproj
index 430ef6b..290272d 100644
--- a/ripser.xcodeproj/project.pbxproj
+++ b/ripser.xcodeproj/project.pbxproj
@@ -208,7 +208,9 @@
 				GCC_PREPROCESSOR_DEFINITIONS = (
 					"$(inherited)",
 					FILE_FORMAT_LOWER_TRIANGULAR_CSV,
+					USE_COEFFICIENTS,
 					STORE_DIAMETERS,
+					PRINT_PERSISTENCE_PAIRS,
 				);
 				PRODUCT_NAME = "$(TARGET_NAME)";
 			};