restructured apparent pairs

2 files changed, 147 insertions, 143 deletions

diff --git a/Makefile b/Makefile
index ab410bd..921e438 100644
--- a/Makefile
+++ b/Makefile
@@ -5,10 +5,10 @@ all: ripser ripser-coeff ripser-debug
 
 
 ripser: ripser.cpp
-	c++ -std=c++11 -Wall ripser.cpp -o ripser -Ofast -D NDEBUG
+	c++ -std=c++11 -Wall ripser.cpp -o ripser -O3 -D NDEBUG
 
 ripser-coeff: ripser.cpp
-	c++ -std=c++11 -Wall ripser.cpp -o ripser-coeff -Ofast -D NDEBUG -D USE_COEFFICIENTS
+	c++ -std=c++11 -Wall ripser.cpp -o ripser-coeff -O3 -D NDEBUG -D USE_COEFFICIENTS
 
 ripser-debug: ripser.cpp
 	c++ -std=c++11 -Wall ripser.cpp -o ripser-debug -g
diff --git a/ripser.cpp b/ripser.cpp
index 1ffa443..6f0fd0a 100644
--- a/ripser.cpp
+++ b/ripser.cpp
@@ -201,13 +201,17 @@ void set_coefficient(diameter_entry_t& p, const coefficient_t c) {
 	set_coefficient(get_entry(p), c);
 }
 
-template <typename Entry> struct greater_diameter_or_smaller_index {
+template <typename Entry> struct greater_diameter_or_smaller_index_comp {
 	bool operator()(const Entry& a, const Entry& b) {
-		return (get_diameter(a) > get_diameter(b)) ||
-		       ((get_diameter(a) == get_diameter(b)) && (get_index(a) < get_index(b)));
+		return greater_diameter_or_smaller_index(a, b);
 	}
 };
 
+template <typename Entry> bool greater_diameter_or_smaller_index(const Entry& a, const Entry& b) {
+	return (get_diameter(a) > get_diameter(b)) ||
+	       ((get_diameter(a) == get_diameter(b)) && (get_index(a) < get_index(b)));
+}
+
 enum compressed_matrix_layout { LOWER_TRIANGULAR, UPPER_TRIANGULAR };
 
 template <compressed_matrix_layout Layout> struct compressed_distance_matrix {
@@ -456,32 +460,31 @@ public:
 	class simplex_boundary_enumerator {
 	private:
 		index_t idx_below, idx_above, j, k;
-        diameter_entry_t simplex;
+		diameter_entry_t simplex;
 		const coefficient_t modulus;
 		const binomial_coeff_table& binomial_coeff;
 		const index_t dim;
 		const ripser& parent;
 
 	public:
-        simplex_boundary_enumerator(const diameter_entry_t _simplex, const index_t _dim,
-                                    const ripser& _parent)
-            : idx_below(get_index(_simplex)), idx_above(0), j(_parent.n - 1), k(_dim),
-              simplex(_simplex), modulus(_parent.modulus), binomial_coeff(_parent.binomial_coeff),
-              dim(_dim), parent(_parent) {}
-
-        simplex_boundary_enumerator(const index_t _dim, const ripser& _parent)
-            : simplex_boundary_enumerator(-1, _dim, _parent) {}
-        
-        void set_simplex(const diameter_entry_t _simplex, const index_t _dim) {
-            idx_below = get_index(_simplex);
-            idx_above = 0;
-            j = parent.n - 1;
-            k = _dim;
-            simplex = _simplex;
-        }
-
-
-        bool has_next() { return (k >= 0); }
+		simplex_boundary_enumerator(const diameter_entry_t _simplex, const index_t _dim,
+		                            const ripser& _parent)
+		    : idx_below(get_index(_simplex)), idx_above(0), j(_parent.n - 1), k(_dim),
+		      simplex(_simplex), modulus(_parent.modulus), binomial_coeff(_parent.binomial_coeff),
+		      dim(_dim), parent(_parent) {}
+
+		simplex_boundary_enumerator(const index_t _dim, const ripser& _parent)
+		    : simplex_boundary_enumerator(-1, _dim, _parent) {}
+
+		void set_simplex(const diameter_entry_t _simplex, const index_t _dim) {
+			idx_below = get_index(_simplex);
+			idx_above = 0;
+			j = parent.n - 1;
+			k = _dim;
+			simplex = _simplex;
+		}
+
+		bool has_next() { return (k >= 0); }
 
 		diameter_entry_t next() {
 			j = parent.get_max_vertex(idx_below, k + 1, j);
@@ -502,57 +505,56 @@ public:
 		}
 	};
 
-	diameter_entry_t get_apparent_facet(const diameter_entry_t simplex, const index_t dim,
-                                        simplex_boundary_enumerator& facets,
-                                        simplex_coboundary_enumerator& cofacets) {
-        facets.set_simplex(simplex, dim);
+	diameter_entry_t get_zero_pivot_facet(const diameter_entry_t simplex, const index_t dim,
+	                                      simplex_boundary_enumerator& facets) {
+		facets.set_simplex(simplex, dim);
 		while (facets.has_next()) {
 			diameter_entry_t facet = facets.next();
-			if (get_diameter(facet) == get_diameter(simplex)) {
-                cofacets.set_simplex(facet, dim - 1);
-				while (cofacets.has_next()) {
-					diameter_entry_t cofacet = cofacets.next();
-					if (get_diameter(cofacet) == get_diameter(simplex))
-						return (get_index(cofacet) == get_index(simplex)) ? facet
-						                                                  : diameter_entry_t(-1);
-				}
-			}
+			if (get_diameter(facet) == get_diameter(simplex)) return facet;
 		}
 		return diameter_entry_t(-1);
 	}
 
-    diameter_entry_t get_apparent_facet(const diameter_entry_t simplex, const index_t dim) {
-        simplex_boundary_enumerator facets(dim, *this);
-        simplex_coboundary_enumerator cofacets(dim - 1, *this);
-        return get_apparent_facet(simplex, dim, facets, cofacets);
-    }
-    
-	diameter_entry_t get_apparent_cofacet(const diameter_entry_t simplex, const index_t dim,
-                                          simplex_coboundary_enumerator& cofacets,
-                                          simplex_boundary_enumerator& facets) {
-        cofacets.set_simplex(simplex, dim);
+	diameter_entry_t get_zero_pivot_cofacet(const diameter_entry_t simplex, const index_t dim,
+	                                        simplex_coboundary_enumerator& cofacets) {
+		cofacets.set_simplex(simplex, dim);
 		while (cofacets.has_next()) {
-            diameter_entry_t cofacet = cofacets.next();
-			if (get_diameter(cofacet) == get_diameter(simplex)) {
-
-                facets.set_simplex(cofacet, dim + 1);
-				while (facets.has_next()) {
-					auto facet = facets.next();
-					if (get_diameter(facet) == get_diameter(simplex))
-						return (get_index(facet) == get_index(simplex)) ? cofacet
-						                                                : diameter_entry_t(-1);
-				}
-			}
+			diameter_entry_t cofacet = cofacets.next();
+			if (get_diameter(cofacet) == get_diameter(simplex)) return cofacet;
 		}
 		return diameter_entry_t(-1);
 	}
 
-    diameter_entry_t get_apparent_cofacet(const diameter_entry_t simplex, const index_t dim) {
-        simplex_coboundary_enumerator cofacets(dim, *this);
-        simplex_boundary_enumerator facets(dim + 1, *this);
-        return get_apparent_cofacet(simplex, dim, cofacets, facets);
-    }
-                
+	diameter_entry_t get_zero_apparent_facet(const diameter_entry_t simplex, const index_t dim,
+	                                         simplex_boundary_enumerator& facets,
+	                                         simplex_coboundary_enumerator& cofacets) {
+		diameter_entry_t facet = get_zero_pivot_facet(simplex, dim, facets);
+		return ((get_index(facet) != -1) &&
+		        (get_index(get_zero_pivot_cofacet(facet, dim - 1, cofacets)) == get_index(simplex)))
+		           ? facet
+		           : diameter_entry_t(-1);
+	}
+
+	diameter_entry_t get_zero_apparent_cofacet(const diameter_entry_t simplex, const index_t dim,
+	                                           simplex_coboundary_enumerator& cofacets,
+	                                           simplex_boundary_enumerator& facets) {
+		diameter_entry_t cofacet = get_zero_pivot_cofacet(simplex, dim, cofacets);
+		return ((get_index(cofacet) != -1) &&
+		        (get_index(get_zero_pivot_facet(cofacet, dim + 1, facets)) == get_index(simplex)))
+		           ? cofacet
+		           : diameter_entry_t(-1);
+	}
+
+	bool is_in_zero_apparent_pair(const diameter_entry_t simplex, const index_t dim,
+	                              simplex_coboundary_enumerator& cofacets_dim,
+	                              simplex_coboundary_enumerator& cofacets_dim_1,
+	                              simplex_boundary_enumerator& facets_dim,
+	                              simplex_boundary_enumerator& facets_dim_1) {
+		return (get_index(get_zero_apparent_cofacet(simplex, dim, cofacets_dim, facets_dim_1)) !=
+		        -1) ||
+		       (get_index(get_zero_apparent_facet(simplex, dim, facets_dim, cofacets_dim_1)) != -1);
+	}
+
 	void assemble_columns_to_reduce(std::vector<diameter_index_t>& simplices,
 	                                std::vector<diameter_index_t>& columns_to_reduce,
 	                                entry_hash_map& pivot_column_index, index_t dim) {
@@ -564,9 +566,10 @@ public:
 
 		columns_to_reduce.clear();
 		std::vector<diameter_index_t> next_simplices;
-        
-        simplex_coboundary_enumerator cofacets(dim - 1, *this), cofacets_dim(dim, *this), cofacets_dim_1(dim - 1, *this);
-        simplex_boundary_enumerator facets_dim(dim, *this), facets_dim_1(dim + 1, *this);
+
+		simplex_coboundary_enumerator cofacets(dim - 1, *this), cofacets_dim(dim, *this),
+		    cofacets_dim_1(dim - 1, *this);
+		simplex_boundary_enumerator facets_dim(dim, *this), facets_dim_1(dim + 1, *this);
 
 		for (diameter_index_t& simplex : simplices) {
 			cofacets.set_simplex(diameter_entry_t(simplex, 1), dim - 1);
@@ -582,12 +585,10 @@ public:
 #endif
 				auto cofacet = cofacets.next();
 				if (get_diameter(cofacet) <= threshold) {
-
 					next_simplices.push_back({get_diameter(cofacet), get_index(cofacet)});
-
-					if ((get_index(get_apparent_cofacet(cofacet, dim, cofacets_dim, facets_dim_1)) == -1) &&
-					    (pivot_column_index.find(get_entry(cofacet)) == pivot_column_index.end()) &&
-					    (get_index(get_apparent_facet(cofacet, dim, facets_dim, cofacets_dim_1)) == -1))
+					if (!is_in_zero_apparent_pair(cofacet, dim, cofacets_dim, cofacets_dim_1,
+					                              facets_dim, facets_dim_1) &&
+					    (pivot_column_index.find(get_entry(cofacet)) == pivot_column_index.end()))
 						columns_to_reduce.push_back({get_diameter(cofacet), get_index(cofacet)});
 				}
 			}
@@ -601,7 +602,7 @@ public:
 #endif
 
 		std::sort(columns_to_reduce.begin(), columns_to_reduce.end(),
-		          greater_diameter_or_smaller_index<diameter_index_t>());
+		          greater_diameter_or_smaller_index<diameter_index_t>);
 #ifdef INDICATE_PROGRESS
 		std::cerr << clear_line << std::flush;
 #endif
@@ -614,13 +615,13 @@ public:
 #endif
 
 		union_find dset(n);
-		
+
 		simplex_coboundary_enumerator cofacets(1, *this);
 		simplex_boundary_enumerator facets(2, *this);
 
 		edges = get_edges();
 		std::sort(edges.rbegin(), edges.rend(),
-		          greater_diameter_or_smaller_index<diameter_index_t>());
+		          greater_diameter_or_smaller_index<diameter_index_t>);
 		std::vector<index_t> vertices_of_edge(2);
 		for (auto e : edges) {
 			get_simplex_vertices(get_index(e), 1, n, vertices_of_edge.rbegin());
@@ -632,7 +633,7 @@ public:
 					std::cout << " [0," << get_diameter(e) << ")" << std::endl;
 #endif
 				dset.link(u, v);
-			} else if (get_index(get_apparent_cofacet(e, 1, cofacets, facets)) == -1)
+			} else if (get_index(get_zero_apparent_cofacet(e, 1, cofacets, facets)) == -1)
 				columns_to_reduce.push_back(e);
 		}
 		std::reverse(columns_to_reduce.begin(), columns_to_reduce.end());
@@ -678,10 +679,9 @@ public:
 	diameter_entry_t init_coboundary_and_get_pivot(const diameter_entry_t simplex,
 	                                               Column& working_coboundary, const index_t& dim,
 	                                               entry_hash_map& pivot_column_index,
-                                                   simplex_coboundary_enumerator& cofacets,
-												   simplex_boundary_enumerator& facets,
-												   simplex_coboundary_enumerator& other_cofacets
-												   ) {
+	                                               simplex_coboundary_enumerator& cofacets,
+	                                               simplex_boundary_enumerator& facets,
+	                                               simplex_coboundary_enumerator& other_cofacets) {
 		bool check_for_emergent_pair = true;
 		cofacet_entries.clear();
 		cofacets.set_simplex(simplex, dim);
@@ -691,7 +691,8 @@ public:
 				cofacet_entries.push_back(cofacet);
 				if (check_for_emergent_pair && (get_diameter(simplex) == get_diameter(cofacet))) {
 					if ((pivot_column_index.find(get_entry(cofacet)) == pivot_column_index.end()) &&
-					    (get_index(get_apparent_facet(cofacet, dim + 1, facets, other_cofacets)) == -1))
+					    (get_index(get_zero_apparent_facet(cofacet, dim + 1, facets,
+					                                       other_cofacets)) == -1))
 						return cofacet;
 					check_for_emergent_pair = false;
 				}
@@ -704,9 +705,9 @@ public:
 	template <typename Column>
 	void add_simplex_coboundary(const diameter_entry_t simplex, const index_t& dim,
 	                            Column& working_reduction_column, Column& working_coboundary,
-                                simplex_coboundary_enumerator& cofacets) {
+	                            simplex_coboundary_enumerator& cofacets) {
 		working_reduction_column.push(simplex);
-        cofacets.set_simplex(simplex, dim);
+		cofacets.set_simplex(simplex, dim);
 		while (cofacets.has_next()) {
 			diameter_entry_t cofacet = cofacets.next();
 			if (get_diameter(cofacet) <= threshold) working_coboundary.push(cofacet);
@@ -718,14 +719,15 @@ public:
 	                    const std::vector<diameter_index_t>& columns_to_reduce,
 	                    const size_t index_column_to_add, const coefficient_t factor,
 	                    const size_t& dim, Column& working_reduction_column,
-	                    Column& working_coboundary,
-                        simplex_coboundary_enumerator& cofacets) {
+	                    Column& working_coboundary, simplex_coboundary_enumerator& cofacets) {
 		diameter_entry_t column_to_add(columns_to_reduce[index_column_to_add], factor);
-		add_simplex_coboundary(column_to_add, dim, working_reduction_column, working_coboundary, cofacets);
+		add_simplex_coboundary(column_to_add, dim, working_reduction_column, working_coboundary,
+		                       cofacets);
 
 		for (diameter_entry_t simplex : reduction_matrix.subrange(index_column_to_add)) {
 			set_coefficient(simplex, get_coefficient(simplex) * factor % modulus);
-			add_simplex_coboundary(simplex, dim, working_reduction_column, working_coboundary, cofacets);
+			add_simplex_coboundary(simplex, dim, working_reduction_column, working_coboundary,
+			                       cofacets);
 		}
 	}
 
@@ -737,10 +739,9 @@ public:
 #endif
 
 		compressed_sparse_matrix<diameter_entry_t> reduction_matrix;
-        
-        simplex_coboundary_enumerator cofacets(dim, *this), other_cofacets(dim, *this);
-        simplex_boundary_enumerator facets(dim + 1, *this);
 
+		simplex_coboundary_enumerator cofacets(dim, *this), other_cofacets(dim, *this);
+		simplex_boundary_enumerator facets(dim + 1, *this);
 
 #ifdef INDICATE_PROGRESS
 		std::chrono::steady_clock::time_point next = std::chrono::steady_clock::now() + time_step;
@@ -754,11 +755,12 @@ public:
 			reduction_matrix.append_column();
 
 			std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>,
-			                    greater_diameter_or_smaller_index<diameter_entry_t>>
+			                    greater_diameter_or_smaller_index_comp<diameter_entry_t>>
 			    working_reduction_column, working_coboundary;
 
 			diameter_entry_t e, pivot = init_coboundary_and_get_pivot(
-			                        column_to_reduce, working_coboundary, dim, pivot_column_index, cofacets, facets, other_cofacets);
+			                        column_to_reduce, working_coboundary, dim, pivot_column_index,
+			                        cofacets, facets, other_cofacets);
 
 			while (true) {
 #ifdef INDICATE_PROGRESS
@@ -780,15 +782,17 @@ public:
 						                  modulus;
 
 						add_coboundary(reduction_matrix, columns_to_reduce, index_column_to_add,
-						               factor, dim, working_reduction_column, working_coboundary, cofacets);
+						               factor, dim, working_reduction_column, working_coboundary,
+						               cofacets);
 
 						pivot = get_pivot(working_coboundary);
-					} else if (get_index(e = get_apparent_facet(pivot, dim + 1, facets, other_cofacets)) != -1) {
+					} else if (get_index(e = get_zero_apparent_facet(pivot, dim + 1, facets,
+					                                                 other_cofacets)) != -1) {
 
 						set_coefficient(e, modulus - get_coefficient(e));
 
-						add_simplex_coboundary(e, dim, working_reduction_column,
-						                       working_coboundary, other_cofacets);
+						add_simplex_coboundary(e, dim, working_reduction_column, working_coboundary,
+						                       other_cofacets);
 
 						pivot = get_pivot(working_coboundary);
 
@@ -851,34 +855,34 @@ public:
 template <> class ripser<compressed_lower_distance_matrix>::simplex_coboundary_enumerator {
 	index_t idx_below, idx_above, j, k;
 	std::vector<index_t> vertices;
-    diameter_entry_t simplex;
+	diameter_entry_t simplex;
 	const coefficient_t modulus;
 	const compressed_lower_distance_matrix& dist;
 	const binomial_coeff_table& binomial_coeff;
-    const ripser& parent;
+	const ripser& parent;
 
 public:
 	simplex_coboundary_enumerator(const diameter_entry_t _simplex, const index_t _dim,
 	                              const ripser& _parent)
 	    : idx_below(get_index(_simplex)), idx_above(0), j(_parent.n - 1), k(_dim + 1),
 	      vertices(_dim + 1), simplex(_simplex), modulus(_parent.modulus), dist(_parent.dist),
-          binomial_coeff(_parent.binomial_coeff), parent(_parent) {
-        if (get_index(_simplex) != -1)
-            parent.get_simplex_vertices(get_index(_simplex), _dim, parent.n, vertices.rbegin());
+	      binomial_coeff(_parent.binomial_coeff), parent(_parent) {
+		if (get_index(_simplex) != -1)
+			parent.get_simplex_vertices(get_index(_simplex), _dim, parent.n, vertices.rbegin());
+	}
+
+	simplex_coboundary_enumerator(const index_t _dim, const ripser& _parent)
+	    : simplex_coboundary_enumerator(-1, _dim, _parent) {}
+
+	void set_simplex(const diameter_entry_t _simplex, const index_t _dim) {
+		idx_below = get_index(_simplex);
+		idx_above = 0;
+		j = parent.n - 1;
+		k = _dim + 1;
+		vertices.resize(_dim + 1);
+		simplex = _simplex;
+		parent.get_simplex_vertices(get_index(_simplex), _dim, parent.n, vertices.rbegin());
 	}
-    
-    simplex_coboundary_enumerator(const index_t _dim, const ripser& _parent)
-        : simplex_coboundary_enumerator(-1, _dim, _parent) {}
-
-    void set_simplex(const diameter_entry_t _simplex, const index_t _dim) {
-        idx_below = get_index(_simplex);
-        idx_above = 0;
-        j = parent.n - 1;
-        k = _dim + 1;
-        vertices.resize(_dim + 1);
-        simplex = _simplex;
-        parent.get_simplex_vertices(get_index(_simplex), _dim, parent.n, vertices.rbegin());
-    }
 
 	bool has_next(bool all_cofacets = true) {
 		return (j >= k && (all_cofacets || binomial_coeff(j, k) > idx_below));
@@ -904,42 +908,42 @@ public:
 template <> class ripser<sparse_distance_matrix>::simplex_coboundary_enumerator {
 	index_t idx_below, idx_above, k;
 	std::vector<index_t> vertices;
-    diameter_entry_t simplex;
+	diameter_entry_t simplex;
 	const coefficient_t modulus;
 	const sparse_distance_matrix& dist;
 	const binomial_coeff_table& binomial_coeff;
 	std::vector<std::vector<index_diameter_t>::const_reverse_iterator> neighbor_it;
 	std::vector<std::vector<index_diameter_t>::const_reverse_iterator> neighbor_end;
 	index_diameter_t neighbor;
-    const ripser& parent;
+	const ripser& parent;
 
 public:
-    simplex_coboundary_enumerator(const diameter_entry_t _simplex, const index_t _dim,
-                                  const ripser& _parent)
-    : idx_below(get_index(_simplex)), idx_above(0), k(_dim + 1), vertices(_dim + 1),
-    simplex(_simplex), modulus(_parent.modulus), dist(_parent.dist),
-    binomial_coeff(_parent.binomial_coeff),
-    neighbor_it(_dim + 1), neighbor_end(_dim + 1), parent(_parent) {
-        if (get_index(_simplex) != -1) set_simplex(_simplex, _dim);
-    }
-    
-    simplex_coboundary_enumerator(const index_t _dim, const ripser& _parent)
-        : simplex_coboundary_enumerator(-1, _dim, _parent) {}
-
-    void set_simplex(const diameter_entry_t _simplex, const index_t _dim) {
-        idx_below = get_index(_simplex);
-        idx_above = 0;
-        k = _dim + 1;
-        simplex = _simplex;
-        
-        parent.get_simplex_vertices(idx_below, _dim, parent.n, vertices.rbegin());
-
-        for (index_t i = 0; i <= _dim; ++i) {
-            auto v = vertices[i];
-            neighbor_it[i] = dist.neighbors[v].rbegin();
-            neighbor_end[i] = dist.neighbors[v].rend();
-        }
-   }
+	simplex_coboundary_enumerator(const diameter_entry_t _simplex, const index_t _dim,
+	                              const ripser& _parent)
+	    : idx_below(get_index(_simplex)), idx_above(0), k(_dim + 1), vertices(_dim + 1),
+	      simplex(_simplex), modulus(_parent.modulus), dist(_parent.dist),
+	      binomial_coeff(_parent.binomial_coeff), neighbor_it(_dim + 1), neighbor_end(_dim + 1),
+	      parent(_parent) {
+		if (get_index(_simplex) != -1) set_simplex(_simplex, _dim);
+	}
+
+	simplex_coboundary_enumerator(const index_t _dim, const ripser& _parent)
+	    : simplex_coboundary_enumerator(-1, _dim, _parent) {}
+
+	void set_simplex(const diameter_entry_t _simplex, const index_t _dim) {
+		idx_below = get_index(_simplex);
+		idx_above = 0;
+		k = _dim + 1;
+		simplex = _simplex;
+
+		parent.get_simplex_vertices(idx_below, _dim, parent.n, vertices.rbegin());
+
+		for (index_t i = 0; i <= _dim; ++i) {
+			auto v = vertices[i];
+			neighbor_it[i] = dist.neighbors[v].rbegin();
+			neighbor_end[i] = dist.neighbors[v].rend();
+		}
+	}
 
 	bool has_next(bool all_cofacets = true) {
 		for (auto &it0 = neighbor_it[0], &end0 = neighbor_end[0]; it0 != end0; ++it0) {