consolidation of branches

1 files changed, 244 insertions, 203 deletions

diff --git a/ripser.cpp b/ripser.cpp
index 8f2d509..9e291ca 100644
--- a/ripser.cpp
+++ b/ripser.cpp
@@ -141,7 +141,8 @@ public:
 	diameter_entry_t(const diameter_index_t& _diameter_index, coefficient_t _coefficient)
 	    : std::pair<value_t, entry_t>(get_diameter(_diameter_index),
 	                                  make_entry(get_index(_diameter_index), _coefficient)) {}
-	diameter_entry_t(const diameter_index_t& _diameter_index) : diameter_entry_t(_diameter_index, 1) {}
+	diameter_entry_t(const diameter_index_t& _diameter_index)
+	    : diameter_entry_t(_diameter_index, 1) {}
 };
 
 const entry_t& get_entry(const diameter_entry_t& p) { return p.second; }
@@ -354,8 +355,8 @@ public:
 	}
 };
 
-class ripser {
-	compressed_lower_distance_matrix dist;
+template <typename DistanceMatrix> class ripser {
+	DistanceMatrix dist;
 	index_t n, dim_max;
 	value_t threshold;
 	float ratio;
@@ -363,9 +364,10 @@ class ripser {
 	const binomial_coeff_table binomial_coeff;
 	std::vector<coefficient_t> multiplicative_inverse;
 	mutable std::vector<index_t> vertices;
+	mutable std::vector<diameter_entry_t> coface_entries;
 
 public:
-	ripser(compressed_lower_distance_matrix&& _dist, index_t _dim_max, value_t _threshold, float _ratio,
+	ripser(DistanceMatrix&& _dist, index_t _dim_max, value_t _threshold, float _ratio,
 	       coefficient_t _modulus)
 	    : dist(std::move(_dist)), n(dist.size()),
 	      dim_max(std::min(_dim_max, index_t(dist.size() - 2))), threshold(_threshold),
@@ -415,83 +417,84 @@ public:
 		return diam;
 	}
 
-	class simplex_coboundary_enumerator {
-	private:
-		index_t idx_below, idx_above, v, k;
-		std::vector<index_t> vertices;
-		const diameter_entry_t simplex;
-		const coefficient_t modulus;
-		const compressed_lower_distance_matrix& dist;
-		const binomial_coeff_table& binomial_coeff;
-
-	public:
-		simplex_coboundary_enumerator(const diameter_entry_t _simplex, index_t _dim,
-		                              const ripser& parent)
-		    : idx_below(get_index(_simplex)), idx_above(0), v(parent.n - 1), k(_dim + 1),
-		      vertices(_dim + 1), simplex(_simplex), modulus(parent.modulus), dist(parent.dist),
-		      binomial_coeff(parent.binomial_coeff) {
-			parent.get_simplex_vertices(get_index(_simplex), _dim, parent.n, vertices.begin());
-		}
+	class simplex_coboundary_enumerator;
 
-		bool has_next() {
-			while ((v != -1) && (binomial_coeff(v, k) <= idx_below)) {
-				idx_below -= binomial_coeff(v, k);
-				idx_above += binomial_coeff(v, k + 1);
-				--v;
-				--k;
-				assert(k != -1);
-			}
-			return v != -1;
-		}
-
-		diameter_entry_t next() {
-			value_t coface_diameter = get_diameter(simplex);
-			for (index_t w : vertices) coface_diameter = std::max(coface_diameter, dist(v, w));
-			index_t coface_index = idx_above + binomial_coeff(v--, k + 1) + idx_below;
-			coefficient_t coface_coefficient =
-			    (k & 1 ? -1 + modulus : 1) * get_coefficient(simplex) % modulus;
-			return diameter_entry_t(coface_diameter, coface_index, coface_coefficient);
-		}
-	};
+	void assemble_columns_to_reduce(std::vector<diameter_index_t>& simplices,
+	                                std::vector<diameter_index_t>& columns_to_reduce,
+	                                hash_map<index_t, index_t>& pivot_column_index, index_t dim);
 
-	void compute_barcodes();
+	void compute_dim_0_pairs(std::vector<diameter_index_t>& edges,
+	                         std::vector<diameter_index_t>& columns_to_reduce) {
+		union_find dset(n);
 
-	void assemble_columns_to_reduce(std::vector<diameter_index_t>& columns_to_reduce,
-	                                hash_map<index_t, index_t>& pivot_column_index, index_t dim) {
-		index_t num_simplices = binomial_coeff(n, dim + 1);
+		edges = get_edges();
 
-		columns_to_reduce.clear();
+		std::sort(edges.rbegin(), edges.rend(),
+		          greater_diameter_or_smaller_index<diameter_index_t>());
 
-#ifdef INDICATE_PROGRESS
-		std::cout << "\033[K"
-		          << "assembling " << num_simplices << " columns" << std::flush << "\r";
+#ifdef PRINT_PERSISTENCE_PAIRS
+		std::cout << "persistence intervals in dim 0:" << std::endl;
 #endif
 
-		for (index_t index = 0; index < num_simplices; ++index) {
-			if (pivot_column_index.find(index) == pivot_column_index.end()) {
-				value_t diameter = compute_diameter(index, dim);
-				if (diameter <= threshold)
-					columns_to_reduce.push_back(std::make_pair(diameter, index));
-#ifdef INDICATE_PROGRESS
-				if ((index + 1) % 1000000 == 0)
-					std::cout << "\033[K"
-					          << "assembled " << columns_to_reduce.size() << " out of "
-					          << (index + 1) << "/" << num_simplices << " columns" << std::flush
-					          << "\r";
+		std::vector<index_t> vertices_of_edge(2);
+		for (auto e : edges) {
+			vertices_of_edge.clear();
+			get_simplex_vertices(get_index(e), 1, n, std::back_inserter(vertices_of_edge));
+			index_t u = dset.find(vertices_of_edge[0]), v = dset.find(vertices_of_edge[1]);
+
+			if (u != v) {
+#ifdef PRINT_PERSISTENCE_PAIRS
+				if (get_diameter(e) != 0)
+					std::cout << " [0," << get_diameter(e) << ")" << std::endl;
 #endif
-			}
+				dset.link(u, v);
+			} else
+				columns_to_reduce.push_back(e);
 		}
+		std::reverse(columns_to_reduce.begin(), columns_to_reduce.end());
 
-#ifdef INDICATE_PROGRESS
-		std::cout << "\033[K"
-		          << "sorting " << num_simplices << " columns" << std::flush << "\r";
+#ifdef PRINT_PERSISTENCE_PAIRS
+		for (index_t i = 0; i < n; ++i)
+			if (dset.find(i) == i) std::cout << " [0, )" << std::endl << std::flush;
 #endif
+	}
 
-		std::sort(columns_to_reduce.begin(), columns_to_reduce.end(),
-		          greater_diameter_or_smaller_index<diameter_index_t>());
-#ifdef INDICATE_PROGRESS
-		std::cout << "\033[K";
+	template <typename Column, typename Iterator>
+	diameter_entry_t add_coboundary_and_get_pivot(Iterator column_begin, Iterator column_end,
+	                                              coefficient_t factor_column_to_add,
+#ifdef ASSEMBLE_REDUCTION_MATRIX
+	                                              Column& working_reduction_column,
 #endif
+	                                              Column& working_coboundary, const index_t& dim,
+	                                              hash_map<index_t, index_t>& pivot_column_index,
+	                                              bool& might_be_apparent_pair) {
+		for (auto it = column_begin; it != column_end; ++it) {
+			diameter_entry_t simplex = *it;
+			set_coefficient(simplex, get_coefficient(simplex) * factor_column_to_add % modulus);
+
+#ifdef ASSEMBLE_REDUCTION_MATRIX
+			working_reduction_column.push(simplex);
+#endif
+
+			coface_entries.clear();
+			simplex_coboundary_enumerator cofaces(simplex, dim, *this);
+			while (cofaces.has_next()) {
+				diameter_entry_t coface = cofaces.next();
+				if (get_diameter(coface) <= threshold) {
+					coface_entries.push_back(coface);
+					if (might_be_apparent_pair && (get_diameter(simplex) == get_diameter(coface))) {
+						if (pivot_column_index.find(get_index(coface)) ==
+						    pivot_column_index.end()) {
+							return coface;
+						}
+						might_be_apparent_pair = false;
+					}
+				}
+			}
+			for (auto coface : coface_entries) working_coboundary.push(coface);
+		}
+
+		return get_pivot(working_coboundary, modulus);
 	}
 
 	void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce,
@@ -511,7 +514,8 @@ public:
 
 		std::vector<diameter_entry_t> coface_entries;
 
-		for (index_t index_column_to_reduce = 0; index_column_to_reduce < columns_to_reduce.size(); ++index_column_to_reduce) {
+		for (index_t index_column_to_reduce = 0; index_column_to_reduce < columns_to_reduce.size();
+		     ++index_column_to_reduce) {
 			auto column_to_reduce = columns_to_reduce[index_column_to_reduce];
 
 #ifdef ASSEMBLE_REDUCTION_MATRIX
@@ -529,14 +533,14 @@ public:
 #ifdef INDICATE_PROGRESS
 			if ((index_column_to_reduce + 1) % 1000000 == 0)
 				std::cout << "\033[K"
-				          << "reducing column " << index_column_to_reduce + 1 << "/" << columns_to_reduce.size()
-				          << " (diameter " << diameter << ")" << std::flush << "\r";
+				          << "reducing column " << index_column_to_reduce + 1 << "/"
+				          << columns_to_reduce.size() << " (diameter " << diameter << ")"
+				          << std::flush << "\r";
 #endif
 
 			index_t index_column_to_add = index_column_to_reduce;
-			
-			diameter_entry_t pivot;
 
+			diameter_entry_t pivot;
 
 			// start with factor 1 in order to initialize working_coboundary
 			// with the coboundary of the simplex with index column_to_reduce
@@ -552,8 +556,7 @@ public:
 
 			bool might_be_apparent_pair = (index_column_to_reduce == index_column_to_add);
 
-			do {
-
+			while (true) {
 #ifdef ASSEMBLE_REDUCTION_MATRIX
 #ifdef USE_COEFFICIENTS
 				auto reduction_column_begin = reduction_matrix.cbegin(index_column_to_add),
@@ -571,39 +574,17 @@ public:
 				auto reduction_column_begin = &reduction_matrix[index_column_to_add],
 				     reduction_column_end = &reduction_matrix[index_column_to_add] + 1;
 #else
-				auto reduction_column_begin = &columns_to_reduce[index_column_to_add], reduction_column_end = &columns_to_reduce[index_column_to_add] + 1;
+				auto reduction_column_begin = &columns_to_reduce[index_column_to_add],
+				     reduction_column_end = &columns_to_reduce[index_column_to_add] + 1;
 #endif
 #endif
 
-				for (auto it = reduction_column_begin; it != reduction_column_end; ++it) {
-					diameter_entry_t simplex = *it;
-					set_coefficient(simplex, get_coefficient(simplex) * factor_column_to_add % modulus);
-
+				pivot = add_coboundary_and_get_pivot(
+				    reduction_column_begin, reduction_column_end, factor_column_to_add,
 #ifdef ASSEMBLE_REDUCTION_MATRIX
-					working_reduction_column.push(simplex);
-#endif
-
-					coface_entries.clear();
-					simplex_coboundary_enumerator cofaces(simplex, dim, *this);
-					while (cofaces.has_next()) {
-						diameter_entry_t coface = cofaces.next();
-						if (get_diameter(coface) <= threshold) {
-							coface_entries.push_back(coface);
-							if (might_be_apparent_pair &&
-							    (get_diameter(simplex) == get_diameter(coface))) {
-								if (pivot_column_index.find(get_index(coface)) ==
-								    pivot_column_index.end()) {
-									pivot = coface;
-									goto found_persistence_pair;
-								}
-								might_be_apparent_pair = false;
-							}
-						}
-					}
-					for (auto coface : coface_entries) working_coboundary.push(coface);
-				}
-
-				pivot = get_pivot(working_coboundary, modulus);
+				    working_reduction_column,
+#endif
+				    working_coboundary, dim, pivot_column_index, might_be_apparent_pair);
 
 				if (get_index(pivot) != -1) {
 					auto pair = pivot_column_index.find(get_index(pivot));
@@ -611,69 +592,174 @@ public:
 					if (pair != pivot_column_index.end()) {
 						index_column_to_add = pair->second;
 						factor_column_to_add = modulus - get_coefficient(pivot);
-						continue;
-					}
-				} else {
+					} else {
 #ifdef PRINT_PERSISTENCE_PAIRS
+						value_t death = get_diameter(pivot);
+						if (diameter != death) {
 #ifdef INDICATE_PROGRESS
-					std::cout << "\033[K";
-#endif
-					std::cout << " [" << diameter << ", )" << std::endl << std::flush;
+							std::cout << "\033[K";
 #endif
-					break;
-				}
-
-			found_persistence_pair:
-#ifdef PRINT_PERSISTENCE_PAIRS
-				value_t death = get_diameter(pivot);
-				if (death > diameter * ratio) {
-#ifdef INDICATE_PROGRESS
-					std::cout << "\033[K";
-#endif
-					std::cout << " [" << diameter << "," << death << ")" << std::endl << std::flush;
-				}
+							std::cout << " [" << diameter << "," << death << ")" << std::endl
+							          << std::flush;
+						}
 #endif
-
-				pivot_column_index.insert(std::make_pair(get_index(pivot), index_column_to_reduce));
+						pivot_column_index.insert(
+						    std::make_pair(get_index(pivot), index_column_to_reduce));
 
 #ifdef USE_COEFFICIENTS
-				const coefficient_t inverse = multiplicative_inverse[get_coefficient(pivot)];
+						const coefficient_t inverse =
+						    multiplicative_inverse[get_coefficient(pivot)];
 #endif
 
 #ifdef ASSEMBLE_REDUCTION_MATRIX
-// replace current column of reduction_matrix (with a single diagonal 1 entry)
-// by reduction_column (possibly with a different entry on the diagonal)
+						// replace current column of reduction_matrix (with a single diagonal 1
+						// entry) by reduction_column (possibly with a different entry on the
+						// diagonal)
 #ifdef USE_COEFFICIENTS
-				reduction_matrix.pop_back();
+						reduction_matrix.pop_back();
 #else
-				pop_pivot(working_reduction_column, modulus);
+						pop_pivot(working_reduction_column, modulus);
 #endif
 
-				while (true) {
-					diameter_entry_t e = pop_pivot(working_reduction_column, modulus);
-					if (get_index(e) == -1) break;
+						while (true) {
+							diameter_entry_t e = pop_pivot(working_reduction_column, modulus);
+							if (get_index(e) == -1) break;
 #ifdef USE_COEFFICIENTS
-					set_coefficient(e, inverse * get_coefficient(e) % modulus);
-					assert(get_coefficient(e) > 0);
+							set_coefficient(e, inverse * get_coefficient(e) % modulus);
+							assert(get_coefficient(e) > 0);
 #endif
-					reduction_matrix.push_back(e);
-				}
+							reduction_matrix.push_back(e);
+						}
 #else
 #ifdef USE_COEFFICIENTS
-				reduction_matrix.pop_back();
-				reduction_matrix.push_back(diameter_entry_t(column_to_reduce, inverse));
+						reduction_matrix.pop_back();
+						reduction_matrix.push_back(diameter_entry_t(column_to_reduce, inverse));
 #endif
 #endif
-				break;
-			} while (true);
+						break;
+					}
+				} else {
+#ifdef PRINT_PERSISTENCE_PAIRS
+					std::cout << " [" << diameter << ", )" << std::endl << std::flush;
+#endif
+					break;
+				}
+			}
 		}
 
 #ifdef INDICATE_PROGRESS
 		std::cout << "\033[K";
 #endif
 	}
+
+	std::vector<diameter_index_t> get_edges();
+
+	void compute_barcodes() {
+
+		std::vector<diameter_index_t> simplices, columns_to_reduce;
+
+		compute_dim_0_pairs(simplices, columns_to_reduce);
+
+		for (index_t dim = 1; dim <= dim_max; ++dim) {
+			hash_map<index_t, index_t> pivot_column_index;
+			pivot_column_index.reserve(columns_to_reduce.size());
+
+			compute_pairs(columns_to_reduce, pivot_column_index, dim);
+
+			if (dim < dim_max) {
+				assemble_columns_to_reduce(simplices, columns_to_reduce, pivot_column_index,
+				                           dim + 1);
+			}
+		}
+	}
+};
+
+template <> class ripser<compressed_lower_distance_matrix>::simplex_coboundary_enumerator {
+private:
+	index_t idx_below, idx_above, v, k;
+	std::vector<index_t> vertices;
+	const diameter_entry_t simplex;
+	const coefficient_t modulus;
+	const compressed_lower_distance_matrix& dist;
+	const binomial_coeff_table& binomial_coeff;
+
+public:
+	simplex_coboundary_enumerator(const diameter_entry_t _simplex, index_t _dim,
+	                              const ripser& parent)
+	    : idx_below(get_index(_simplex)), idx_above(0), v(parent.n - 1), k(_dim + 1),
+	      vertices(_dim + 1), simplex(_simplex), modulus(parent.modulus), dist(parent.dist),
+	      binomial_coeff(parent.binomial_coeff) {
+		parent.get_simplex_vertices(get_index(_simplex), _dim, parent.n, vertices.begin());
+	}
+
+	bool has_next() {
+		while ((v != -1) && (binomial_coeff(v, k) <= idx_below)) {
+			idx_below -= binomial_coeff(v, k);
+			idx_above += binomial_coeff(v, k + 1);
+			--v;
+			--k;
+			assert(k != -1);
+		}
+		return v != -1;
+	}
+
+	diameter_entry_t next() {
+		value_t coface_diameter = get_diameter(simplex);
+		for (index_t w : vertices) coface_diameter = std::max(coface_diameter, dist(v, w));
+		index_t coface_index = idx_above + binomial_coeff(v--, k + 1) + idx_below;
+		coefficient_t coface_coefficient =
+		    (k & 1 ? -1 + modulus : 1) * get_coefficient(simplex) % modulus;
+		return diameter_entry_t(coface_diameter, coface_index, coface_coefficient);
+	}
 };
 
+template <> std::vector<diameter_index_t> ripser<compressed_lower_distance_matrix>::get_edges() {
+	std::vector<diameter_index_t> edges;
+	for (index_t index = binomial_coeff(n, 2); index-- > 0;) {
+		value_t diameter = compute_diameter(index, 1);
+		if (diameter <= threshold) edges.push_back(std::make_pair(diameter, index));
+	}
+	return edges;
+}
+
+template <>
+void ripser<compressed_lower_distance_matrix>::assemble_columns_to_reduce(
+    std::vector<diameter_index_t>& simplices, std::vector<diameter_index_t>& columns_to_reduce,
+    hash_map<index_t, index_t>& pivot_column_index, index_t dim) {
+	index_t num_simplices = binomial_coeff(n, dim + 1);
+
+	columns_to_reduce.clear();
+
+#ifdef INDICATE_PROGRESS
+	std::cout << "\033[K"
+	          << "assembling " << num_simplices << " columns" << std::flush << "\r";
+#endif
+
+	for (index_t index = 0; index < num_simplices; ++index) {
+		if (pivot_column_index.find(index) == pivot_column_index.end()) {
+			value_t diameter = compute_diameter(index, dim);
+			if (diameter <= threshold) columns_to_reduce.push_back(std::make_pair(diameter, index));
+#ifdef INDICATE_PROGRESS
+			if ((index + 1) % 1000000 == 0)
+				std::cout << "\033[K"
+				          << "assembled " << columns_to_reduce.size() << " out of " << (index + 1)
+				          << "/" << num_simplices << " columns" << std::flush << "\r";
+#endif
+		}
+	}
+
+#ifdef INDICATE_PROGRESS
+	std::cout << "\033[K"
+	          << "sorting " << num_simplices << " columns" << std::flush << "\r";
+#endif
+
+	std::sort(columns_to_reduce.begin(), columns_to_reduce.end(),
+	          greater_diameter_or_smaller_index<diameter_index_t>());
+#ifdef INDICATE_PROGRESS
+	std::cout << "\033[K";
+#endif
+}
+
 enum file_format {
 	LOWER_DISTANCE_MATRIX,
 	UPPER_DISTANCE_MATRIX,
@@ -908,74 +994,29 @@ int main(int argc, char** argv) {
 
 	compressed_lower_distance_matrix dist = read_file(filename ? file_stream : std::cin, format);
 
-	std::cout << "distance matrix with " << dist.size() << " points" << std::endl;
+	value_t min = std::numeric_limits<value_t>::infinity(),
+	        max = -std::numeric_limits<value_t>::infinity(), max_finite = max;
+	int num_edges = 0;
 
-	value_t min = std::numeric_limits<value_t>::infinity(), max = -std::numeric_limits<value_t>::infinity();
-	
-	for (auto d: dist.distances) {
-		if (d != std::numeric_limits<value_t>::infinity() ) {
-			min = std::min(min, d);
-			max = std::max(max, d);
-		} else {
-			threshold = std::min(threshold, std::numeric_limits<value_t>::max());
-		}
+	value_t enclosing_radius = std::numeric_limits<value_t>::infinity();
+	for (index_t i = 0; i < dist.size(); ++i) {
+		value_t r_i = -std::numeric_limits<value_t>::infinity();
+		for (index_t j = 0; j < dist.size(); ++j) r_i = std::max(r_i, dist(i, j));
+		enclosing_radius = std::min(enclosing_radius, r_i);
 	}
-	
-	std::cout << "value range: [" << min << "," << max << "]"
-	          << std::endl;
-
-	ripser(std::move(dist), dim_max, threshold, ratio,  modulus).compute_barcodes();
-}
-
-void ripser::compute_barcodes() {
-
-	std::vector<diameter_index_t> columns_to_reduce;
-
-	{
-		union_find dset(n);
-		std::vector<diameter_index_t> edges;
-		for (index_t index = binomial_coeff(n, 2); index-- > 0;) {
-			value_t diameter = compute_diameter(index, 1);
-			if (diameter <= threshold) edges.push_back(std::make_pair(diameter, index));
-		}
-		std::sort(edges.rbegin(), edges.rend(),
-		          greater_diameter_or_smaller_index<diameter_index_t>());
-
-#ifdef PRINT_PERSISTENCE_PAIRS
-		std::cout << "persistence intervals in dim 0:" << std::endl;
-#endif
 
-		std::vector<index_t> vertices_of_edge(2);
-		for (auto e : edges) {
-			vertices_of_edge.clear();
-			get_simplex_vertices(get_index(e), 1, n, std::back_inserter(vertices_of_edge));
-			index_t u = dset.find(vertices_of_edge[0]), v = dset.find(vertices_of_edge[1]);
+	if (threshold == std::numeric_limits<value_t>::max()) threshold = enclosing_radius;
 
-			if (u != v) {
-#ifdef PRINT_PERSISTENCE_PAIRS
-				if (get_diameter(e) != 0)
-					std::cout << " [0," << get_diameter(e) << ")" << std::endl;
-#endif
-				dset.link(u, v);
-			} else
-				columns_to_reduce.push_back(e);
-		}
-		std::reverse(columns_to_reduce.begin(), columns_to_reduce.end());
-
-#ifdef PRINT_PERSISTENCE_PAIRS
-		for (index_t i = 0; i < n; ++i)
-			if (dset.find(i) == i) std::cout << " [0, )" << std::endl << std::flush;
-#endif
+	for (auto d : dist.distances) {
+		min = std::min(min, d);
+		max = std::max(max, d);
+		max_finite = d != std::numeric_limits<value_t>::infinity() ? std::max(max, d) : max_finite;
+		if (d <= threshold) ++num_edges;
 	}
 
-	for (index_t dim = 1; dim <= dim_max; ++dim) {
-		hash_map<index_t, index_t> pivot_column_index;
-		pivot_column_index.reserve(columns_to_reduce.size());
+	std::cout << "value range: [" << min << "," << max_finite << "]" << std::endl;
 
-		compute_pairs(columns_to_reduce, pivot_column_index, dim);
-
-		if (dim < dim_max) {
-			assemble_columns_to_reduce(columns_to_reduce, pivot_column_index, dim + 1);
-		}
-	}
+	std::cout << "distance matrix with " << dist.size() << " points" << std::endl;
+	ripser<compressed_lower_distance_matrix>(std::move(dist), dim_max, threshold, ratio, modulus)
+	    .compute_barcodes();
 }