code reformatting

2 files changed, 259 insertions, 212 deletions

diff --git a/.clang-format b/.clang-format
index 598d7c2..1975b19 100644
--- a/.clang-format
+++ b/.clang-format
@@ -1,7 +1,7 @@
 BasedOnStyle: LLVM
 IndentWidth: 4
 TabWidth: 4
-ColumnLimit: 120
+ColumnLimit: 100
 AccessModifierOffset: -4
 AllowShortIfStatementsOnASingleLine: true
 AllowShortLoopsOnASingleLine: true
diff --git a/ripser.cpp b/ripser.cpp
index 052c07d..02bfac9 100644
--- a/ripser.cpp
+++ b/ripser.cpp
@@ -52,6 +52,7 @@ typedef int16_t coefficient_t;
 
 class binomial_coeff_table {
 	std::vector<std::vector<index_t>> B;
+
 public:
 	binomial_coeff_table(index_t n, index_t k) : B(n + 1) {
 		for (index_t i = 0; i <= n; i++) {
@@ -87,19 +88,21 @@ std::vector<coefficient_t> multiplicative_inverse_vector(const coefficient_t m)
 	return inverse;
 }
 
-
 #ifdef USE_COEFFICIENTS
 struct __attribute__((packed)) entry_t {
 	index_t index : 8 * (sizeof(index_t) - sizeof(coefficient_t));
 	coefficient_t coefficient;
-	entry_t(index_t _index, coefficient_t _coefficient) : index(_index), coefficient(_coefficient) {}
+	entry_t(index_t _index, coefficient_t _coefficient)
+	    : index(_index), coefficient(_coefficient) {}
 	entry_t(index_t _index) : index(_index), coefficient(1) {}
 	entry_t() : index(0), coefficient(1) {}
 };
 
 static_assert(sizeof(entry_t) == sizeof(index_t), "size of entry_t is not the same as index_t");
 
-entry_t make_entry(index_t _index, coefficient_t _coefficient) { return entry_t(_index, _coefficient); }
+entry_t make_entry(index_t _index, coefficient_t _coefficient) {
+	return entry_t(_index, _coefficient);
+}
 index_t get_index(entry_t e) { return e.index; }
 index_t get_coefficient(entry_t e) { return e.coefficient; }
 void set_coefficient(entry_t& e, const coefficient_t c) { e.coefficient = c; }
@@ -149,9 +152,13 @@ public:
 const entry_t& get_entry(const diameter_entry_t& p) { return p.second; }
 entry_t& get_entry(diameter_entry_t& p) { return p.second; }
 const index_t get_index(const diameter_entry_t& p) { return get_index(get_entry(p)); }
-const coefficient_t get_coefficient(const diameter_entry_t& p) { return get_coefficient(get_entry(p)); }
+const coefficient_t get_coefficient(const diameter_entry_t& p) {
+	return get_coefficient(get_entry(p));
+}
 const value_t& get_diameter(const diameter_entry_t& p) { return p.first; }
-void set_coefficient(diameter_entry_t& p, const coefficient_t c) { set_coefficient(get_entry(p), c); }
+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 {
 	bool operator()(const Entry& a, const Entry& b) {
@@ -205,11 +212,13 @@ template <> void compressed_distance_matrix<UPPER_TRIANGULAR>::init_rows() {
 	}
 }
 
-template <> value_t compressed_distance_matrix<UPPER_TRIANGULAR>::operator()(index_t i, index_t j) const {
+template <>
+value_t compressed_distance_matrix<UPPER_TRIANGULAR>::operator()(index_t i, index_t j) const {
 	return i == j ? 0 : i > j ? rows[j][i] : rows[i][j];
 }
 
-template <> value_t compressed_distance_matrix<LOWER_TRIANGULAR>::operator()(index_t i, index_t j) const {
+template <>
+value_t compressed_distance_matrix<LOWER_TRIANGULAR>::operator()(index_t i, index_t j) const {
 	return i == j ? 0 : i < j ? rows[j][i] : rows[i][j];
 }
 
@@ -220,12 +229,13 @@ class euclidean_distance_matrix {
 public:
 	std::vector<std::vector<value_t>> points;
 
-	euclidean_distance_matrix(std::vector<std::vector<value_t>>&& _points) : points(std::move(_points)) {}
+	euclidean_distance_matrix(std::vector<std::vector<value_t>>&& _points)
+	    : points(std::move(_points)) {}
 
 	value_t operator()(const index_t i, const index_t j) const {
-		return std::sqrt(std::inner_product(points[i].begin(), points[i].end(), points[j].begin(), value_t(),
-		                                    std::plus<value_t>(),
-		                                    [](value_t u, value_t v) { return (u - v) * (u - v); }));
+		return std::sqrt(std::inner_product(
+		    points[i].begin(), points[i].end(), points[j].begin(), value_t(), std::plus<value_t>(),
+		    [](value_t u, value_t v) { return (u - v) * (u - v); }));
 	}
 
 	size_t size() const { return points.size(); }
@@ -350,7 +360,8 @@ public:
 	}
 };
 
-template <typename Heap> void push_entry(Heap& column, index_t i, coefficient_t c, value_t diameter) {
+template <typename Heap>
+void push_entry(Heap& column, index_t i, coefficient_t c, value_t diameter) {
 	entry_t e = make_entry(i, c);
 	column.push(std::make_pair(diameter, e));
 }
@@ -365,11 +376,12 @@ class ripser {
 	mutable std::vector<index_t> vertices;
 
 public:
-	ripser(compressed_lower_distance_matrix&& _dist, index_t _dim_max, value_t _threshold, coefficient_t _modulus)
-	    : dist(std::move(_dist)), dim_max(std::min(_dim_max, index_t(dist.size() - 2))), n(dist.size()), threshold(_threshold),
-	      modulus(_modulus), binomial_coeff(n, dim_max + 2),
+	ripser(compressed_lower_distance_matrix&& _dist, index_t _dim_max, value_t _threshold,
+	       coefficient_t _modulus)
+	    : dist(std::move(_dist)), dim_max(std::min(_dim_max, index_t(dist.size() - 2))),
+	      n(dist.size()), threshold(_threshold), modulus(_modulus), binomial_coeff(n, dim_max + 2),
 	      multiplicative_inverse(multiplicative_inverse_vector(_modulus)) {}
-	
+
 	index_t get_next_vertex(index_t& v, const index_t idx, const index_t k) const {
 		if (binomial_coeff(v, k) > idx) {
 			index_t count = v;
@@ -387,10 +399,10 @@ public:
 		assert(binomial_coeff(v, k) <= idx && binomial_coeff(v + 1, k) > idx);
 		return v;
 	}
-	
+
 	template <typename OutputIterator>
 	OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index_t v,
-										OutputIterator out) const {
+	                                    OutputIterator out) const {
 		--v;
 		for (index_t k = dim + 1; k > 0; --k) {
 			get_next_vertex(v, idx, k);
@@ -399,7 +411,7 @@ public:
 		}
 		return out;
 	}
-	
+
 	value_t compute_diameter(const index_t index, index_t dim) const {
 		value_t diam = -std::numeric_limits<value_t>::infinity();
 
@@ -407,10 +419,12 @@ public:
 		get_simplex_vertices(index, dim, dist.size(), std::back_inserter(vertices));
 
 		for (index_t i = 0; i <= dim; ++i)
-			for (index_t j = 0; j < i; ++j) { diam = std::max(diam, dist(vertices[i], vertices[j])); }
+			for (index_t j = 0; j < i; ++j) {
+				diam = std::max(diam, dist(vertices[i], vertices[j]));
+			}
 		return diam;
 	}
-	
+
 	class simplex_coboundary_enumerator {
 	private:
 		index_t idx_below, idx_above, v, k;
@@ -419,14 +433,16 @@ public:
 		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) {
+		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);
@@ -437,276 +453,248 @@ public:
 			}
 			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;
+			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>& columns_to_reduce,
-	                                hash_map<index_t, index_t>& pivot_column_index, index_t dim) {
-		index_t num_simplices = binomial_coeff(n, dim + 1);
+	                                hash_map<index_t, index_t>& pivot_column_index, index_t dim);
+
+	void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce,
+	                   hash_map<index_t, index_t>& pivot_column_index, index_t dim);
+
+	void compute_barcodes();
+};
 
-		columns_to_reduce.clear();
+void ripser::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);
+
+	columns_to_reduce.clear();
 
 #ifdef INDICATE_PROGRESS
-		std::cout << "\033[K"
-		          << "assembling " << num_simplices << " columns" << std::flush << "\r";
+	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));
+	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) % 1000 == 0)
-					std::cout << "\033[K"
-					          << "assembled " << columns_to_reduce.size() << " out of " << (index + 1) << "/"
-					          << num_simplices << " columns" << std::flush << "\r";
+			if ((index + 1) % 1000 == 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";
+	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>());
+	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";
+	std::cout << "\033[K";
 #endif
-	}
+}
 
-	void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce, hash_map<index_t, index_t>& pivot_column_index,
-	                   index_t dim) {
+void ripser::compute_pairs(std::vector<diameter_index_t>& columns_to_reduce,
+                           hash_map<index_t, index_t>& pivot_column_index, index_t dim) {
 
 #ifdef PRINT_PERSISTENCE_PAIRS
-		std::cout << "persistence intervals in dim " << dim << ":" << std::endl;
+	std::cout << "persistence intervals in dim " << dim << ":" << std::endl;
 #endif
 
 #ifdef ASSEMBLE_REDUCTION_MATRIX
-		compressed_sparse_matrix<diameter_entry_t> reduction_coefficients;
+	compressed_sparse_matrix<diameter_entry_t> reduction_coefficients;
 #else
 #ifdef USE_COEFFICIENTS
-		std::vector<diameter_entry_t> reduction_coefficients;
+	std::vector<diameter_entry_t> reduction_coefficients;
 #endif
 #endif
 
-		std::vector<diameter_entry_t> coface_entries;
+	std::vector<diameter_entry_t> coface_entries;
 
-		for (index_t i = 0; i < columns_to_reduce.size(); ++i) {
-			auto column_to_reduce = columns_to_reduce[i];
+	for (index_t i = 0; i < columns_to_reduce.size(); ++i) {
+		auto column_to_reduce = columns_to_reduce[i];
 
 #ifdef ASSEMBLE_REDUCTION_MATRIX
-			std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>,
-			                    greater_diameter_or_smaller_index<diameter_entry_t>>
-			    reduction_column;
+		std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>,
+		                    greater_diameter_or_smaller_index<diameter_entry_t>>
+		    reduction_column;
 #endif
 
-			std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>,
-			                    greater_diameter_or_smaller_index<diameter_entry_t>>
-			    working_coboundary;
+		std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>,
+		                    greater_diameter_or_smaller_index<diameter_entry_t>>
+		    working_coboundary;
 
-			value_t diameter = get_diameter(column_to_reduce);
+		value_t diameter = get_diameter(column_to_reduce);
 
 #ifdef INDICATE_PROGRESS
-			if ((i + 1) % 1000 == 0)
-				std::cout << "\033[K"
-				          << "reducing column " << i + 1 << "/" << columns_to_reduce.size() << " (diameter " << diameter
-				          << ")" << std::flush << "\r";
+		if ((i + 1) % 1000 == 0)
+			std::cout << "\033[K"
+			          << "reducing column " << i + 1 << "/" << columns_to_reduce.size()
+			          << " (diameter " << diameter << ")" << std::flush << "\r";
 #endif
 
-			index_t j = i;
+		index_t j = i;
 
-			// 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
-			diameter_entry_t pivot(0, -1, -1 + modulus);
+		// 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
+		diameter_entry_t pivot(0, -1, -1 + modulus);
 
 #ifdef ASSEMBLE_REDUCTION_MATRIX
-			// initialize reduction_coefficients as identity matrix
-			reduction_coefficients.append_column();
+		// initialize reduction_coefficients as identity matrix
+		reduction_coefficients.append_column();
 #endif
 #ifdef USE_COEFFICIENTS
-			reduction_coefficients.push_back(diameter_entry_t(column_to_reduce, 1));
+		reduction_coefficients.push_back(diameter_entry_t(column_to_reduce, 1));
 #endif
 
-			bool might_be_apparent_pair = (i == j);
+		bool might_be_apparent_pair = (i == j);
 
-			do {
-				const coefficient_t factor = modulus - get_coefficient(pivot);
+		do {
+			const coefficient_t factor = modulus - get_coefficient(pivot);
 
 #ifdef ASSEMBLE_REDUCTION_MATRIX
 #ifdef USE_COEFFICIENTS
-				auto coeffs_begin = reduction_coefficients.cbegin(j), coeffs_end = reduction_coefficients.cend(j);
+			auto coeffs_begin = reduction_coefficients.cbegin(j),
+			     coeffs_end = reduction_coefficients.cend(j);
 #else
-				std::vector<diameter_entry_t> coeffs;
-				coeffs.push_back(columns_to_reduce[j]);
-				for (auto it = reduction_coefficients.cbegin(j); it != reduction_coefficients.cend(j); ++it)
-					coeffs.push_back(*it);
-				auto coeffs_begin = coeffs.begin(), coeffs_end = coeffs.end();
+			std::vector<diameter_entry_t> coeffs;
+			coeffs.push_back(columns_to_reduce[j]);
+			for (auto it = reduction_coefficients.cbegin(j); it != reduction_coefficients.cend(j);
+			     ++it)
+				coeffs.push_back(*it);
+			auto coeffs_begin = coeffs.begin(), coeffs_end = coeffs.end();
 #endif
 #else
 #ifdef USE_COEFFICIENTS
-				auto coeffs_begin = &reduction_coefficients[j], coeffs_end = &reduction_coefficients[j] + 1;
+			auto coeffs_begin = &reduction_coefficients[j],
+			     coeffs_end = &reduction_coefficients[j] + 1;
 #else
-				auto coeffs_begin = &columns_to_reduce[j], coeffs_end = &columns_to_reduce[j] + 1;
+			auto coeffs_begin = &columns_to_reduce[j], coeffs_end = &columns_to_reduce[j] + 1;
 #endif
 #endif
 
-				for (auto it = coeffs_begin; it != coeffs_end; ++it) {
-					diameter_entry_t simplex = *it;
-					set_coefficient(simplex, get_coefficient(simplex) * factor % modulus);
+			for (auto it = coeffs_begin; it != coeffs_end; ++it) {
+				diameter_entry_t simplex = *it;
+				set_coefficient(simplex, get_coefficient(simplex) * factor % modulus);
 
 #ifdef ASSEMBLE_REDUCTION_MATRIX
-					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;
+				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);
 				}
+				for (auto coface : coface_entries) working_coboundary.push(coface);
+			}
 
-				pivot = get_pivot(working_coboundary, modulus);
+			pivot = get_pivot(working_coboundary, modulus);
 
-				if (get_index(pivot) != -1) {
-					auto pair = pivot_column_index.find(get_index(pivot));
+			if (get_index(pivot) != -1) {
+				auto pair = pivot_column_index.find(get_index(pivot));
 
-					if (pair != pivot_column_index.end()) {
-						j = pair->second;
-						continue;
-					}
-				} else {
+				if (pair != pivot_column_index.end()) {
+					j = pair->second;
+					continue;
+				}
+			} else {
 #ifdef PRINT_PERSISTENCE_PAIRS
 #ifdef INDICATE_PROGRESS
-					std::cout << "\033[K";
+				std::cout << "\033[K";
 #endif
-					std::cout << " [" << diameter << ", )" << std::endl << std::flush;
+				std::cout << " [" << diameter << ", )" << std::endl << std::flush;
 #endif
-					break;
-				}
+				break;
+			}
 
-			found_persistence_pair:
+		found_persistence_pair:
 #ifdef PRINT_PERSISTENCE_PAIRS
-				value_t death = get_diameter(pivot);
-				if (diameter != death) {
+			value_t death = get_diameter(pivot);
+			if (diameter != death) {
 #ifdef INDICATE_PROGRESS
-					std::cout << "\033[K";
+				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), i));
+			pivot_column_index.insert(std::make_pair(get_index(pivot), i));
 
 #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_coefficients (with a single diagonal 1 entry)
 // by reduction_column (possibly with a different entry on the diagonal)
 #ifdef USE_COEFFICIENTS
-				reduction_coefficients.pop_back();
+			reduction_coefficients.pop_back();
 #else
-				pop_pivot(reduction_column, modulus);
+			pop_pivot(reduction_column, modulus);
 #endif
 
-				while (true) {
-					diameter_entry_t e = pop_pivot(reduction_column, modulus);
-					if (get_index(e) == -1) break;
+			while (true) {
+				diameter_entry_t e = pop_pivot(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_coefficients.push_back(e);
-				}
+				reduction_coefficients.push_back(e);
+			}
 #else
 #ifdef USE_COEFFICIENTS
-				reduction_coefficients.pop_back();
-				reduction_coefficients.push_back(diameter_entry_t(column_to_reduce, inverse));
+			reduction_coefficients.pop_back();
+			reduction_coefficients.push_back(diameter_entry_t(column_to_reduce, inverse));
 #endif
 #endif
-				break;
-			} while (true);
-		}
-
-#ifdef INDICATE_PROGRESS
-		std::cout << "\033[K";
-#endif
+			break;
+		} while (true);
 	}
 
-	void 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 (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;
+#ifdef INDICATE_PROGRESS
+	std::cout << "\033[K";
 #endif
-		}
-
-		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(columns_to_reduce, pivot_column_index, dim + 1); }
-		}
-	}
+enum file_format {
+	LOWER_DISTANCE_MATRIX,
+	UPPER_DISTANCE_MATRIX,
+	DISTANCE_MATRIX,
+	POINT_CLOUD,
+	DIPHA,
+	RIPSER
 };
 
-enum file_format { LOWER_DISTANCE_MATRIX, UPPER_DISTANCE_MATRIX, DISTANCE_MATRIX, POINT_CLOUD, DIPHA, RIPSER };
-
 template <typename T> T read(std::istream& s) {
 	T result;
 	s.read(reinterpret_cast<char*>(&result), sizeof(T));
@@ -733,7 +721,8 @@ compressed_lower_distance_matrix read_point_cloud(std::istream& input_stream) {
 
 	index_t n = eucl_dist.size();
 
-	std::cout << "point cloud with " << n << " points in dimension " << eucl_dist.points.front().size() << std::endl;
+	std::cout << "point cloud with " << n << " points in dimension "
+	          << eucl_dist.points.front().size() << std::endl;
 
 	std::vector<value_t> distances;
 
@@ -830,26 +819,30 @@ compressed_lower_distance_matrix read_file(std::istream& input_stream, file_form
 }
 
 void print_usage_and_exit(int exit_code) {
-	std::cerr << "Usage: "
-	          << "ripser "
-	          << "[options] [filename]" << std::endl
-	          << std::endl
-	          << "Options:" << std::endl
-	          << std::endl
-	          << "  --help           print this screen" << std::endl
-	          << "  --format         use the specified file format for the input. Options are:" << std::endl
-	          << "                     lower-distance (lower triangular distance matrix; default)" << std::endl
-	          << "                     upper-distance (upper triangular distance matrix)" << std::endl
-	          << "                     distance       (full distance matrix)" << std::endl
-	          << "                     point-cloud    (point cloud in Euclidean space)" << std::endl
-	          << "                     dipha          (distance matrix in DIPHA file format)" << std::endl
-	          << "                     ripser         (distance matrix in Ripser binary file format)" << std::endl
-	          << "  --dim <k>        compute persistent homology up to dimension <k>" << std::endl
-	          << "  --threshold <t>  compute Rips complexes up to diameter <t>" << std::endl
+	std::cerr
+	    << "Usage: "
+	    << "ripser "
+	    << "[options] [filename]" << std::endl
+	    << std::endl
+	    << "Options:" << std::endl
+	    << std::endl
+	    << "  --help           print this screen" << std::endl
+	    << "  --format         use the specified file format for the input. Options are:"
+	    << std::endl
+	    << "                     lower-distance (lower triangular distance matrix; default)"
+	    << std::endl
+	    << "                     upper-distance (upper triangular distance matrix)" << std::endl
+	    << "                     distance       (full distance matrix)" << std::endl
+	    << "                     point-cloud    (point cloud in Euclidean space)" << std::endl
+	    << "                     dipha          (distance matrix in DIPHA file format)" << std::endl
+	    << "                     ripser         (distance matrix in Ripser binary file format)"
+	    << std::endl
+	    << "  --dim <k>        compute persistent homology up to dimension <k>" << std::endl
+	    << "  --threshold <t>  compute Rips complexes up to diameter <t>" << std::endl
 #ifdef USE_COEFFICIENTS
-	          << "  --modulus <p>    compute homology with coefficients in the prime field Z/<p>Z"
+	    << "  --modulus <p>    compute homology with coefficients in the prime field Z/<p>Z"
 #endif
-	          << std::endl;
+	    << std::endl;
 
 	exit(exit_code);
 }
@@ -923,7 +916,61 @@ int main(int argc, char** argv) {
 	std::cout << "distance matrix with " << dist.size() << " points" << std::endl;
 
 	auto value_range = std::minmax_element(dist.distances.begin(), dist.distances.end());
-	std::cout << "value range: [" << *value_range.first << "," << *value_range.second << "]" << std::endl;
+	std::cout << "value range: [" << *value_range.first << "," << *value_range.second << "]"
+	          << std::endl;
 
 	ripser(std::move(dist), dim_max, threshold, 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 (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 (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(columns_to_reduce, pivot_column_index, dim + 1);
+		}
+	}
+}