Merge tag 'v1.2.1' into debian/sid

5 files changed, 253 insertions, 95 deletions

diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 0000000..b85c2fb
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,3 @@
+[submodule "robin-hood-hashing"]
+	path = robin-hood-hashing
+	url = https://github.com/martinus/robin-hood-hashing.git
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/README.md b/README.md
index 51029bf..7472c33 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@
 # Ripser
 
-Copyright © 2015–2019 [Ulrich Bauer].
+Copyright © 2015–2021 [Ulrich Bauer].
 
 
 ### Description
@@ -39,7 +39,7 @@ Ripser's efficiency is based on a few important concepts and principles, buildin
 
 
 ### Version
-[Latest release][latest-release]: 1.1 (August 2019)
+[Latest release][latest-release]: 1.2.1 (March 2021)
 
 
 ### Building
@@ -61,12 +61,12 @@ Ripser supports several compile-time options. They are switched on by defining t
   - `USE_COEFFICIENTS`: enable support for coefficients in a prime field
   - `INDICATE_PROGRESS`: indicate the current progress in the console
   - `PRINT_PERSISTENCE_PAIRS`: output the computed persistence pairs (enabled by default in the code; comment out to disable)
-  - `USE_GOOGLE_HASHMAP`: enable support for Google's [sparsehash] data structure; may further reduce memory footprint
+  - `USE_ROBINHOOD_HASHMAP`: enable support for Martin Ankerl's [robinhoodhash] data structure; may further reduce memory footprint
 
-For example, to build Ripser with support for Google's hashmap:
+For example, to build Ripser with support for Martin Ankerl's robin hood hashmap:
 
 ```sh
-$ c++ -std=c++11 ripser.cpp -o ripser -Ofast -D NDEBUG -D USE_GOOGLE_HASHMAP
+$ c++ -std=c++11 ripser.cpp -o ripser -O3 -D NDEBUG -D USE_ROBINHOOD_HASHMAP
 ```
 
 A Makefile is provided with some variants of the above options. Use `make all` to build them. The default `make` builds a binary with the default options.
@@ -74,12 +74,12 @@ A Makefile is provided with some variants of the above options. Use `make all` t
 The input is given either in a file whose name is passed as an argument, or through stdin. The following options are supported at the command line:
 
   - `--format`: use the specified file format for the input.  The following formats are supported:
-    - `lower-distance` (default if no format is specified): lower triangular distance matrix; a comma (or whitespace, or other non-numerical character) separated list of the distance matrix entries below the diagonal, sorted lexicographically by row index, then column index.
+    - `lower-distance`: lower triangular distance matrix; a comma (or whitespace, or other non-numerical character) separated list of the distance matrix entries below the diagonal, sorted lexicographically by row index, then column index.
     - `upper-distance`: upper triangular distance matrix; similar to the previous, but for the entries above the diagonal; suitable for output from the MATLAB functions `pdist` or  `seqpdist`, exported to a CSV file.
-    - `distance`: full distance matrix; similar to the above, but for all entries of the distance matrix. One line per row of the matrix; only the part below the diagonal is actually read.
+    - `distance`  (default if no format is specified): full distance matrix; similar to the above, but for all entries of the distance matrix. One line per row of the matrix; only the part below the diagonal is actually read.
     - `dipha`: DIPHA distance matrix as described on the [DIPHA] website.
     - `point-cloud`: point cloud; a comma (or whitespace, or other non-numerical character) separated list of coordinates of the points in some Euclidean space, one point per line.
-    - `binary`: lower distance matrix in binary file format; a sequence of the distance matrix entries below the diagonal in 32 bit double format (IEEE 754, little endian).
+    - `binary`: lower distance matrix in binary file format; a sequence of the distance matrix entries below the diagonal in 32 bit float format (IEEE 754, single, little endian).
     - `sparse`: sparse triplet format; a whitespace separated list of entries of a sparse distance matrix, one entry per line, each of the form *i j d(i,j)* specifying the distance between points *i* and *j*.  Each pair of points should appear in the file at most once.
   - `--dim k`: compute persistent homology up to dimension *k*.
   - `--threshold t`: compute Rips complexes up to diameter *t*.
@@ -105,9 +105,9 @@ If you use Ripser in your research or if you want to give a reference to Ripser
 @misc{1908.02518,
 	Author = {Ulrich Bauer},
 	Title = {Ripser: efficient computation of Vietoris-Rips persistence barcodes},
-	Month = Aug,
-	Year = {2019},
-	Eprint = {1908.02518},
+	Month = Feb,
+	Year = {2021},
+	Eprint = {1908.02518v2},
 	Note = {Preprint}
 }
 ```
@@ -127,5 +127,5 @@ Ripser is licensed under the [MIT] license (`COPYING.txt`), with an extra clause
 [PHAT]: <http://git.io/dipha>
 [Perseus]: <http://www.sas.upenn.edu/~vnanda/perseus/>
 [GUDHI]: <http://gudhi.gforge.inria.fr>
-[sparsehash]: <https://github.com/sparsehash/sparsehash>
-[MIT]: <https://opensource.org/licenses/mit-license.php>
\ No newline at end of file
+[robinhoodhash]: <https://github.com/martinus/robin-hood-hashing>
+[MIT]: <https://opensource.org/licenses/mit-license.php>
diff --git a/ripser.cpp b/ripser.cpp
index cd546e8..cfa01a8 100644
--- a/ripser.cpp
+++ b/ripser.cpp
@@ -4,7 +4,7 @@
 
  MIT License
 
- Copyright (c) 2015–2019 Ulrich Bauer
+ Copyright (c) 2015–2021 Ulrich Bauer
 
  Permission is hereby granted, free of charge, to any person obtaining a copy
  of this software and associated documentation files (the "Software"), to deal
@@ -41,7 +41,7 @@
 //#define INDICATE_PROGRESS
 #define PRINT_PERSISTENCE_PAIRS
 
-//#define USE_GOOGLE_HASHMAP
+//#define USE_ROBINHOOD_HASHMAP
 
 #include <algorithm>
 #include <cassert>
@@ -54,17 +54,19 @@
 #include <sstream>
 #include <unordered_map>
 
-#ifdef USE_GOOGLE_HASHMAP
-#include <sparsehash/dense_hash_map>
+#ifdef USE_ROBINHOOD_HASHMAP
+
+#include "robin-hood-hashing/src/include/robin_hood.h"
+
 template <class Key, class T, class H, class E>
-class hash_map : public google::dense_hash_map<Key, T, H, E> {
-public:
-	explicit hash_map() : google::dense_hash_map<Key, T, H, E>() { this->set_empty_key(-1); }
-	inline void reserve(size_t hint) { this->resize(hint); }
-};
+using hash_map = robin_hood::unordered_map<Key, T, H, E>;
+template <class Key> using hash = robin_hood::hash<Key>;
+
 #else
-template <class Key, class T, class H, class E>
-class hash_map : public std::unordered_map<Key, T, H, E> {};
+
+template <class Key, class T, class H, class E> using hash_map = std::unordered_map<Key, T, H, E>;
+template <class Key> using hash = std::hash<Key>;
+
 #endif
 
 typedef float value_t;
@@ -182,6 +184,9 @@ struct diameter_entry_t : std::pair<value_t, entry_t> {
 	diameter_entry_t(const diameter_index_t& _diameter_index, coefficient_t _coefficient)
 	    : diameter_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(get_diameter(_diameter_index),
+	                       make_entry(get_index(_diameter_index), 0)) {}
 	diameter_entry_t(const index_t& _index) : diameter_entry_t(0, _index, 0) {}
 };
 
@@ -196,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 {
@@ -263,9 +272,6 @@ struct sparse_distance_matrix {
 
 	index_t num_edges;
 
-	mutable std::vector<std::vector<index_diameter_t>::const_reverse_iterator> neighbor_it;
-	mutable std::vector<std::vector<index_diameter_t>::const_reverse_iterator> neighbor_end;
-
 	sparse_distance_matrix(std::vector<std::vector<index_diameter_t>>&& _neighbors,
 	                       index_t _num_edges)
 	    : neighbors(std::move(_neighbors)), num_edges(_num_edges) {}
@@ -276,12 +282,23 @@ struct sparse_distance_matrix {
 
 		for (size_t i = 0; i < size(); ++i)
 			for (size_t j = 0; j < size(); ++j)
-				if (i != j && mat(i, j) <= threshold) {
-					++num_edges;
-					neighbors[i].push_back({j, mat(i, j)});
+				if (i != j) {
+					auto d = mat(i, j);
+					if (d <= threshold) {
+						++num_edges;
+						neighbors[i].push_back({j, d});
+					}
 				}
 	}
 
+	value_t operator()(const index_t i, const index_t j) const {
+		auto neighbor =
+		    std::lower_bound(neighbors[i].begin(), neighbors[i].end(), index_diameter_t{j, 0});
+		return (neighbor != neighbors[i].end() && get_index(*neighbor) == j)
+		           ? get_diameter(*neighbor)
+		           : std::numeric_limits<value_t>::infinity();
+	}
+
 	size_t size() const { return neighbors.size(); }
 };
 
@@ -386,11 +403,10 @@ template <typename DistanceMatrix> class ripser {
 	const binomial_coeff_table binomial_coeff;
 	const std::vector<coefficient_t> multiplicative_inverse;
 	mutable std::vector<diameter_entry_t> cofacet_entries;
+	mutable std::vector<index_t> vertices;
 
 	struct entry_hash {
-		std::size_t operator()(const entry_t& e) const {
-			return std::hash<index_t>()(::get_index(e));
-		}
+		std::size_t operator()(const entry_t& e) const { return hash<index_t>()(::get_index(e)); }
 	};
 
 	struct equal_index {
@@ -429,8 +445,111 @@ public:
 		return out;
 	}
 
+	value_t compute_diameter(const index_t index, const index_t dim) const {
+		value_t diam = -std::numeric_limits<value_t>::infinity();
+
+		vertices.resize(dim + 1);
+		get_simplex_vertices(index, dim, dist.size(), vertices.rbegin());
+
+		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]));
+			}
+		return diam;
+	}
+
 	class simplex_coboundary_enumerator;
 
+	class simplex_boundary_enumerator {
+	private:
+		index_t idx_below, idx_above, j, k;
+		diameter_entry_t simplex;
+		index_t dim;
+		const coefficient_t modulus;
+		const binomial_coeff_table& binomial_coeff;
+		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),
+		      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;
+			dim = _dim;
+		}
+
+		bool has_next() { return (k >= 0); }
+
+		diameter_entry_t next() {
+			j = parent.get_max_vertex(idx_below, k + 1, j);
+
+			index_t face_index = idx_above - binomial_coeff(j, k + 1) + idx_below;
+
+			value_t face_diameter = parent.compute_diameter(face_index, dim - 1);
+
+			coefficient_t face_coefficient =
+			    (k & 1 ? -1 + modulus : 1) * get_coefficient(simplex) % modulus;
+
+			idx_below -= binomial_coeff(j, k + 1);
+			idx_above += binomial_coeff(j, k);
+
+			--k;
+
+			return diameter_entry_t(face_diameter, face_index, face_coefficient);
+		}
+	};
+
+	diameter_entry_t get_zero_pivot_facet(const diameter_entry_t simplex, const index_t dim) {
+		static simplex_boundary_enumerator facets(0, *this);
+		facets.set_simplex(simplex, dim);
+		while (facets.has_next()) {
+			diameter_entry_t facet = facets.next();
+			if (get_diameter(facet) == get_diameter(simplex)) return facet;
+		}
+		return diameter_entry_t(-1);
+	}
+
+	diameter_entry_t get_zero_pivot_cofacet(const diameter_entry_t simplex, const index_t dim) {
+		static simplex_coboundary_enumerator cofacets(*this);
+		cofacets.set_simplex(simplex, dim);
+		while (cofacets.has_next()) {
+			diameter_entry_t cofacet = cofacets.next();
+			if (get_diameter(cofacet) == get_diameter(simplex)) return cofacet;
+		}
+		return diameter_entry_t(-1);
+	}
+
+	diameter_entry_t get_zero_apparent_facet(const diameter_entry_t simplex, const index_t dim) {
+		diameter_entry_t facet = get_zero_pivot_facet(simplex, dim);
+		return ((get_index(facet) != -1) &&
+		        (get_index(get_zero_pivot_cofacet(facet, dim - 1)) == get_index(simplex)))
+		           ? facet
+		           : diameter_entry_t(-1);
+	}
+
+	diameter_entry_t get_zero_apparent_cofacet(const diameter_entry_t simplex, const index_t dim) {
+		diameter_entry_t cofacet = get_zero_pivot_cofacet(simplex, dim);
+		return ((get_index(cofacet) != -1) &&
+		        (get_index(get_zero_pivot_facet(cofacet, dim + 1)) == get_index(simplex)))
+		           ? cofacet
+		           : diameter_entry_t(-1);
+	}
+
+	bool is_in_zero_apparent_pair(const diameter_entry_t simplex, const index_t dim) {
+		return (get_index(get_zero_apparent_cofacet(simplex, dim)) != -1) ||
+		       (get_index(get_zero_apparent_facet(simplex, dim)) != -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) {
@@ -440,12 +559,13 @@ public:
 		std::chrono::steady_clock::time_point next = std::chrono::steady_clock::now() + time_step;
 #endif
 
-		--dim;
 		columns_to_reduce.clear();
 		std::vector<diameter_index_t> next_simplices;
 
+		simplex_coboundary_enumerator cofacets(*this);
+
 		for (diameter_index_t& simplex : simplices) {
-			simplex_coboundary_enumerator cofacets(diameter_entry_t(simplex, 1), dim, *this);
+			cofacets.set_simplex(diameter_entry_t(simplex, 1), dim - 1);
 
 			while (cofacets.has_next(false)) {
 #ifdef INDICATE_PROGRESS
@@ -458,10 +578,9 @@ public:
 #endif
 				auto cofacet = cofacets.next();
 				if (get_diameter(cofacet) <= threshold) {
-
 					next_simplices.push_back({get_diameter(cofacet), get_index(cofacet)});
-
-					if (pivot_column_index.find(get_entry(cofacet)) == pivot_column_index.end())
+					if (!is_in_zero_apparent_pair(cofacet, dim) &&
+					    (pivot_column_index.find(get_entry(cofacet)) == pivot_column_index.end()))
 						columns_to_reduce.push_back({get_diameter(cofacet), get_index(cofacet)});
 				}
 			}
@@ -475,7 +594,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
@@ -491,7 +610,7 @@ public:
 
 		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());
@@ -503,7 +622,7 @@ public:
 					std::cout << " [0," << get_diameter(e) << ")" << std::endl;
 #endif
 				dset.link(u, v);
-			} else
+			} else if (get_index(get_zero_apparent_cofacet(e, 1)) == -1)
 				columns_to_reduce.push_back(e);
 		}
 		std::reverse(columns_to_reduce.begin(), columns_to_reduce.end());
@@ -549,15 +668,17 @@ 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) {
+		static simplex_coboundary_enumerator cofacets(*this);
 		bool check_for_emergent_pair = true;
 		cofacet_entries.clear();
-		simplex_coboundary_enumerator cofacets(simplex, dim, *this);
+		cofacets.set_simplex(simplex, dim);
 		while (cofacets.has_next()) {
 			diameter_entry_t cofacet = cofacets.next();
 			if (get_diameter(cofacet) <= threshold) {
 				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())
+					if ((pivot_column_index.find(get_entry(cofacet)) == pivot_column_index.end()) &&
+					    (get_index(get_zero_apparent_facet(cofacet, dim + 1)) == -1))
 						return cofacet;
 					check_for_emergent_pair = false;
 				}
@@ -570,8 +691,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) {
+		static simplex_coboundary_enumerator cofacets(*this);
 		working_reduction_column.push(simplex);
-		simplex_coboundary_enumerator cofacets(simplex, dim, *this);
+		cofacets.set_simplex(simplex, dim);
 		while (cofacets.has_next()) {
 			diameter_entry_t cofacet = cofacets.next();
 			if (get_diameter(cofacet) <= threshold) working_coboundary.push(cofacet);
@@ -581,8 +703,9 @@ public:
 	template <typename Column>
 	void add_coboundary(compressed_sparse_matrix<diameter_entry_t>& reduction_matrix,
 	                    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) {
+	                    const size_t index_column_to_add, const coefficient_t factor,
+	                    const size_t& dim, Column& working_reduction_column,
+	                    Column& working_coboundary) {
 		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);
 
@@ -600,7 +723,7 @@ public:
 #endif
 
 		compressed_sparse_matrix<diameter_entry_t> reduction_matrix;
-
+		
 #ifdef INDICATE_PROGRESS
 		std::chrono::steady_clock::time_point next = std::chrono::steady_clock::now() + time_step;
 #endif
@@ -613,11 +736,11 @@ 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 pivot = init_coboundary_and_get_pivot(
-			    column_to_reduce, working_coboundary, dim, pivot_column_index);
+			diameter_entry_t e, pivot = init_coboundary_and_get_pivot(
+			                        column_to_reduce, working_coboundary, dim, pivot_column_index);
 
 			while (true) {
 #ifdef INDICATE_PROGRESS
@@ -642,6 +765,12 @@ public:
 						               factor, dim, working_reduction_column, working_coboundary);
 
 						pivot = get_pivot(working_coboundary);
+					} else if (get_index(e = get_zero_apparent_facet(pivot, dim + 1)) != -1) {
+						set_coefficient(e, modulus - get_coefficient(e));
+
+						add_simplex_coboundary(e, dim, working_reduction_column, working_coboundary);
+
+						pivot = get_pivot(working_coboundary);
 					} else {
 #ifdef PRINT_PERSISTENCE_PAIRS
 						value_t death = get_diameter(pivot);
@@ -699,37 +828,51 @@ public:
 };
 
 template <> class ripser<compressed_lower_distance_matrix>::simplex_coboundary_enumerator {
-	index_t idx_below, idx_above, v, k;
+	index_t idx_below, idx_above, j, k;
 	std::vector<index_t> vertices;
-	const 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;
 
 public:
 	simplex_coboundary_enumerator(const diameter_entry_t _simplex, const 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) {
+	                              const ripser& _parent)
+	    : 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());
+	}
+
+	simplex_coboundary_enumerator(const ripser& _parent) : modulus(_parent.modulus), dist(_parent.dist),
+	binomial_coeff(_parent.binomial_coeff), parent(_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;
+		simplex = _simplex;
+		vertices.resize(_dim + 1);
 		parent.get_simplex_vertices(get_index(_simplex), _dim, parent.n, vertices.rbegin());
 	}
 
 	bool has_next(bool all_cofacets = true) {
-		return (v >= k && (all_cofacets || binomial_coeff(v, k) > idx_below));
+		return (j >= k && (all_cofacets || binomial_coeff(j, k) > idx_below));
 	}
 
 	diameter_entry_t next() {
-		while ((binomial_coeff(v, k) <= idx_below)) {
-			idx_below -= binomial_coeff(v, k);
-			idx_above += binomial_coeff(v, k + 1);
-			--v;
+		while ((binomial_coeff(j, k) <= idx_below)) {
+			idx_below -= binomial_coeff(j, k);
+			idx_above += binomial_coeff(j, k + 1);
+			--j;
 			--k;
 			assert(k != -1);
 		}
 		value_t cofacet_diameter = get_diameter(simplex);
-		for (index_t w : vertices) cofacet_diameter = std::max(cofacet_diameter, dist(v, w));
-		index_t cofacet_index = idx_above + binomial_coeff(v--, k + 1) + idx_below;
+		for (index_t i : vertices) cofacet_diameter = std::max(cofacet_diameter, dist(j, i));
+		index_t cofacet_index = idx_above + binomial_coeff(j--, k + 1) + idx_below;
 		coefficient_t cofacet_coefficient =
 		    (k & 1 ? modulus - 1 : 1) * get_coefficient(simplex) % modulus;
 		return diameter_entry_t(cofacet_diameter, cofacet_index, cofacet_coefficient);
@@ -737,31 +880,43 @@ public:
 };
 
 template <> class ripser<sparse_distance_matrix>::simplex_coboundary_enumerator {
-	const ripser& parent;
 	index_t idx_below, idx_above, k;
 	std::vector<index_t> vertices;
-	const 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;
+	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;
 
 public:
 	simplex_coboundary_enumerator(const diameter_entry_t _simplex, const index_t _dim,
 	                              const ripser& _parent)
-	    : parent(_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(dist.neighbor_it),
-	      neighbor_end(dist.neighbor_end) {
-		neighbor_it.clear();
-		neighbor_end.clear();
+	    : modulus(_parent.modulus), dist(_parent.dist),
+	      binomial_coeff(_parent.binomial_coeff), parent(_parent) {
+		if (get_index(_simplex) != -1) set_simplex(_simplex, _dim);
+	}
 
+	simplex_coboundary_enumerator(const ripser& _parent)
+	    : modulus(_parent.modulus), dist(_parent.dist),
+	binomial_coeff(_parent.binomial_coeff), parent(_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;
+		vertices.resize(_dim + 1);
 		parent.get_simplex_vertices(idx_below, _dim, parent.n, vertices.rbegin());
-		for (auto v : vertices) {
-			neighbor_it.push_back(dist.neighbors[v].rbegin());
-			neighbor_end.push_back(dist.neighbors[v].rend());
+
+		neighbor_it.resize(_dim + 1);
+		neighbor_end.resize(_dim + 1);
+		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();
 		}
 	}
 
@@ -841,7 +996,7 @@ template <typename T> T read(std::istream& input_stream) {
 	return result;
 }
 
-compressed_lower_distance_matrix read_point_cloud(std::istream& input_stream) {
+euclidean_distance_matrix read_point_cloud(std::istream& input_stream) {
 	std::vector<std::vector<value_t>> points;
 
 	std::string line;
@@ -862,11 +1017,7 @@ compressed_lower_distance_matrix read_point_cloud(std::istream& input_stream) {
 	std::cout << "point cloud with " << n << " points in dimension "
 	          << eucl_dist.points.front().size() << std::endl;
 
-	std::vector<value_t> distances;
-	for (int i = 0; i < n; ++i)
-		for (int j = 0; j < i; ++j) distances.push_back(eucl_dist(i, j));
-
-	return compressed_lower_distance_matrix(std::move(distances));
+	return eucl_dist;
 }
 
 sparse_distance_matrix read_sparse_distance_matrix(std::istream& input_stream) {
@@ -1088,6 +1239,10 @@ int main(int argc, char** argv) {
 
 		ripser<sparse_distance_matrix>(std::move(dist), dim_max, threshold, ratio, modulus)
 		    .compute_barcodes();
+	} else if (format == POINT_CLOUD && threshold < std::numeric_limits<value_t>::max()) {
+		sparse_distance_matrix dist(read_point_cloud(filename ? file_stream : std::cin), threshold);
+		ripser<sparse_distance_matrix>(std::move(dist), dim_max, threshold, ratio, modulus)
+				.compute_barcodes();
 	} else {
 		compressed_lower_distance_matrix dist =
 		    read_file(filename ? file_stream : std::cin, format);
@@ -1096,29 +1251,29 @@ int main(int argc, char** argv) {
 		        max = -std::numeric_limits<value_t>::infinity(), max_finite = max;
 		int num_edges = 0;
 
+		value_t enclosing_radius = std::numeric_limits<value_t>::infinity();
 		if (threshold == std::numeric_limits<value_t>::max()) {
-			value_t enclosing_radius = std::numeric_limits<value_t>::infinity();
 			for (size_t i = 0; i < dist.size(); ++i) {
 				value_t r_i = -std::numeric_limits<value_t>::infinity();
 				for (size_t j = 0; j < dist.size(); ++j) r_i = std::max(r_i, dist(i, j));
 				enclosing_radius = std::min(enclosing_radius, r_i);
 			}
-			threshold = enclosing_radius;
 		}
 
 		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 != std::numeric_limits<value_t>::infinity()) max_finite = std::max(max_finite, d);
 			if (d <= threshold) ++num_edges;
 		}
 		std::cout << "value range: [" << min << "," << max_finite << "]" << std::endl;
 
-		if (threshold >= max) {
-			std::cout << "distance matrix with " << dist.size() << " points" << std::endl;
-			ripser<compressed_lower_distance_matrix>(std::move(dist), dim_max, threshold, ratio,
-			                                         modulus)
+		if (threshold == std::numeric_limits<value_t>::max()) {
+			std::cout << "distance matrix with " << dist.size()
+			          << " points, using threshold at enclosing radius " << enclosing_radius
+			          << std::endl;
+			ripser<compressed_lower_distance_matrix>(std::move(dist), dim_max, enclosing_radius,
+			                                         ratio, modulus)
 			    .compute_barcodes();
 		} else {
 			std::cout << "sparse distance matrix with " << dist.size() << " points and "
diff --git a/robin-hood-hashing b/robin-hood-hashing
new file mode 160000
+Subproject 7b00237399cec62ee5795f9013cf5df8de596e6