updated Google hash map support

1 files changed, 108 insertions, 137 deletions

diff --git a/ripser.cpp b/ripser.cpp
index e7c3d0e..7178df8 100644
--- a/ripser.cpp
+++ b/ripser.cpp
@@ -1,6 +1,3 @@
-#define USE_BINARY_SEARCH
-//#define USE_EXPONENTIAL_SEARCH
-
 //#define ASSEMBLE_REDUCTION_MATRIX
 //#define USE_COEFFICIENTS
 
@@ -56,7 +53,7 @@ public:
 		}
 	}
 
-	inline index_t operator()(index_t n, index_t k) const {
+	index_t operator()(index_t n, index_t k) const {
 		assert(n <= n_max);
 		assert(k <= k_max);
 		return B[n][k];
@@ -66,34 +63,22 @@ public:
 std::vector<coefficient_t> multiplicative_inverse_vector(const coefficient_t m) {
 	std::vector<coefficient_t> inverse(m);
 	inverse[1] = 1;
-	for (coefficient_t a = 2; a < m; ++a) {
-		// m = a * (m / a) + m % a
-		// Multipying with inverse(a) * inverse(m % a):
-		// 0 = (m / a) * inverse(m % a) + inverse(a) = 0   (mod m)
-		inverse[a] = m - ((m / a) * inverse[m % a]) % m;
-	}
+	// m = a * (m / a) + m % a
+	// Multipying with inverse(a) * inverse(m % a):
+	// 0 = (m / a) * inverse(m % a) + inverse(a)  (mod m)
+	for (coefficient_t a = 2; a < m; ++a) inverse[a] = m - ((m / a) * inverse[m % a]) % m;
 	return inverse;
 }
 
 template <typename OutputIterator>
-inline OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index_t n,
-                                           const binomial_coeff_table& binomial_coeff,
-                                           OutputIterator out) {
+OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index_t n,
+                                    const binomial_coeff_table& binomial_coeff,
+                                    OutputIterator out) {
 	--n;
 
 	for (index_t k = dim + 1; k > 0; --k) {
-
-#ifdef USE_BINARY_SEARCH
 		if (binomial_coeff(n, k) > idx) {
-			index_t count;
-
-#ifdef USE_EXPONENTIAL_SEARCH
-			for (count = 1; (binomial_coeff(n - count, k) > idx); count = std::min(count << 1, n))
-				;
-#else
-			count = n;
-#endif
-
+			index_t count = n;
 			while (count > 0) {
 				index_t i = n;
 				index_t step = count >> 1;
@@ -105,15 +90,10 @@ inline OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index
 					count = step;
 			}
 		}
-#else
-		while (binomial_coeff(n, k) > idx) --n;
-#endif
-
 		assert(binomial_coeff(n, k) <= idx);
 		assert(binomial_coeff(n + 1, k) > idx);
 
 		*out++ = n;
-
 		idx -= binomial_coeff(n, k);
 	}
 
@@ -138,12 +118,12 @@ struct entry_t {
 	entry_t() : index(0), coefficient(1) {}
 };
 
-inline entry_t make_entry(index_t _index, coefficient_t _coefficient) {
+entry_t make_entry(index_t _index, coefficient_t _coefficient) {
 	return entry_t(_index, _coefficient);
 }
-inline index_t get_index(entry_t e) { return e.index; }
-inline index_t get_coefficient(entry_t e) { return e.coefficient; }
-inline void set_coefficient(entry_t& e, const coefficient_t c) { e.coefficient = c; }
+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; }
 
 bool operator==(const entry_t& e1, const entry_t& e2) {
 	return get_index(e1) == get_index(e2) && get_coefficient(e1) == get_coefficient(e2);
@@ -157,22 +137,22 @@ std::ostream& operator<<(std::ostream& stream, const entry_t& e) {
 #else
 
 typedef index_t entry_t;
-inline const index_t get_index(entry_t i) { return i; }
-inline index_t get_coefficient(entry_t i) { return 1; }
-inline entry_t make_entry(index_t _index, coefficient_t _value) { return entry_t(_index); }
-inline void set_coefficient(index_t& e, const coefficient_t c) { e = c; }
+const index_t get_index(entry_t i) { return i; }
+index_t get_coefficient(entry_t i) { return 1; }
+entry_t make_entry(index_t _index, coefficient_t _value) { return entry_t(_index); }
+void set_coefficient(index_t& e, const coefficient_t c) { e = c; }
 
 #endif
 
-inline const entry_t& get_entry(const entry_t& e) { return e; }
+const entry_t& get_entry(const entry_t& e) { return e; }
 
 template <typename Entry> struct smaller_index {
 	bool operator()(const Entry& a, const Entry& b) { return get_index(a) < get_index(b); }
 };
 
 typedef std::pair<value_t, index_t> diameter_index_t;
-inline value_t get_diameter(diameter_index_t i) { return i.first; }
-inline index_t get_index(diameter_index_t i) { return i.second; }
+value_t get_diameter(diameter_index_t i) { return i.first; }
+index_t get_index(diameter_index_t i) { return i.second; }
 
 class diameter_entry_t : public std::pair<value_t, entry_t> {
 public:
@@ -181,28 +161,27 @@ public:
 	diameter_entry_t() : diameter_entry_t(0) {}
 };
 
-inline const entry_t& get_entry(const diameter_entry_t& p) { return p.second; }
-inline entry_t& get_entry(diameter_entry_t& p) { return p.second; }
-inline const index_t get_index(const diameter_entry_t& p) { return get_index(get_entry(p)); }
-inline const coefficient_t get_coefficient(const diameter_entry_t& p) {
+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));
 }
-inline const value_t& get_diameter(const diameter_entry_t& p) { return p.first; }
-inline void set_coefficient(diameter_entry_t& p, const coefficient_t c) {
+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);
 }
-inline diameter_entry_t make_diameter_entry(value_t _diameter, index_t _index,
-                                            coefficient_t _coefficient) {
+diameter_entry_t make_diameter_entry(value_t _diameter, index_t _index,
+                                     coefficient_t _coefficient) {
 	return std::make_pair(_diameter, make_entry(_index, _coefficient));
 }
-inline diameter_entry_t make_diameter_entry(diameter_index_t _diameter_index,
-                                            coefficient_t _coefficient) {
+diameter_entry_t make_diameter_entry(diameter_index_t _diameter_index, coefficient_t _coefficient) {
 	return std::make_pair(get_diameter(_diameter_index),
 	                      make_entry(get_index(_diameter_index), _coefficient));
 }
 
 template <typename Entry> struct greater_diameter_or_smaller_index {
-	inline bool operator()(const Entry& a, const Entry& b) {
+	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)));
 	}
@@ -222,7 +201,7 @@ public:
 	                           const binomial_coeff_table& _binomial_coeff)
 	    : dist(_dist), dim(_dim), vertices(_dim + 1), binomial_coeff(_binomial_coeff){};
 
-	inline value_t diameter(const index_t index) const {
+	value_t diameter(const index_t index) const {
 		value_t diam = 0;
 		get_simplex_vertices(index, dim, dist.size(), binomial_coeff, vertices.begin());
 
@@ -233,7 +212,7 @@ public:
 		return diam;
 	}
 
-	inline bool operator()(const index_t a, const index_t b) const {
+	bool operator()(const index_t a, const index_t b) const {
 		assert(a < binomial_coeff(dist.size(), dim + 1));
 		assert(b < binomial_coeff(dist.size(), dim + 1));
 
@@ -249,7 +228,7 @@ public:
 		return (a_diam == b_diam) && (a < b);
 	}
 
-	template <typename Entry> inline bool operator()(const Entry& a, const Entry& b) const {
+	template <typename Entry> bool operator()(const Entry& a, const Entry& b) const {
 		return operator()(get_index(a), get_index(b));
 	}
 };
@@ -260,12 +239,12 @@ private:
 	const binomial_coeff_table& binomial_coeff;
 
 public:
-	inline simplex_coboundary_enumerator(index_t _idx, index_t _dim, index_t _n,
-	                                     const binomial_coeff_table& _binomial_coeff)
+	simplex_coboundary_enumerator(index_t _idx, index_t _dim, index_t _n,
+	                              const binomial_coeff_table& _binomial_coeff)
 	    : idx(_idx), modified_idx(_idx), dim(_dim), k(dim + 1), v(_n - 1),
 	      binomial_coeff(_binomial_coeff) {}
 
-	inline bool has_next() {
+	bool has_next() {
 		while ((v != -1) && (binomial_coeff(v, k) <= idx)) {
 			idx -= binomial_coeff(v, k);
 			modified_idx += binomial_coeff(v, k + 1) - binomial_coeff(v, k);
@@ -276,7 +255,7 @@ public:
 		return v != -1;
 	}
 
-	inline std::pair<entry_t, index_t> next() {
+	std::pair<entry_t, index_t> next() {
 		auto result =
 		    std::make_pair(make_entry(modified_idx + binomial_coeff(v, k + 1), k & 1 ? -1 : 1), v);
 		--v;
@@ -288,9 +267,9 @@ class distance_matrix {
 public:
 	typedef value_t value_type;
 	std::vector<std::vector<value_t>> distances;
-	inline value_t operator()(const index_t a, const index_t b) const { return distances[a][b]; }
+	value_t operator()(const index_t a, const index_t b) const { return distances[a][b]; }
 
-	inline size_t size() const { return distances.size(); }
+	size_t size() const { return distances.size(); }
 };
 
 enum compressed_matrix_layout { LOWER_TRIANGULAR, UPPER_TRIANGULAR };
@@ -322,9 +301,9 @@ public:
 		}
 	}
 
-	inline value_t operator()(const index_t i, const index_t j) const;
+	value_t operator()(const index_t i, const index_t j) const;
 
-	inline size_t size() const { return rows.size(); }
+	size_t size() const { return rows.size(); }
 };
 
 template <> void compressed_distance_matrix_adapter<LOWER_TRIANGULAR>::init_rows() {
@@ -344,8 +323,8 @@ template <> void compressed_distance_matrix_adapter<UPPER_TRIANGULAR>::init_rows
 }
 
 template <>
-inline value_t compressed_distance_matrix_adapter<UPPER_TRIANGULAR>::
-operator()(const index_t i, const index_t j) const {
+value_t compressed_distance_matrix_adapter<UPPER_TRIANGULAR>::operator()(const index_t i,
+                                                                         const index_t j) const {
 	if (i < j)
 		return rows[i][j];
 	else if (i > j)
@@ -355,8 +334,8 @@ operator()(const index_t i, const index_t j) const {
 }
 
 template <>
-inline value_t compressed_distance_matrix_adapter<LOWER_TRIANGULAR>::
-operator()(const index_t i, const index_t j) const {
+value_t compressed_distance_matrix_adapter<LOWER_TRIANGULAR>::operator()(const index_t i,
+                                                                         const index_t j) const {
 	if (i > j)
 		return rows[i][j];
 	else if (i < j)
@@ -370,7 +349,7 @@ typedef compressed_distance_matrix_adapter<LOWER_TRIANGULAR>
 typedef compressed_distance_matrix_adapter<UPPER_TRIANGULAR>
     compressed_upper_distance_matrix_adapter;
 
-template <typename Heap> inline diameter_entry_t pop_pivot(Heap& column, coefficient_t modulus) {
+template <typename Heap> diameter_entry_t pop_pivot(Heap& column, coefficient_t modulus) {
 	if (column.empty())
 		return diameter_entry_t(-1);
 	else {
@@ -407,9 +386,9 @@ template <typename Heap> inline diameter_entry_t pop_pivot(Heap& column, coeffic
 	}
 }
 
-template <typename Heap> inline diameter_entry_t get_pivot(Heap& column, coefficient_t modulus) {
+template <typename Heap> diameter_entry_t get_pivot(Heap& column, coefficient_t modulus) {
 	diameter_entry_t result = pop_pivot(column, modulus);
-	if (get_index(result) != -1) { column.push(result); }
+	if (get_index(result) != -1) column.push(result);
 	return result;
 }
 
@@ -438,44 +417,44 @@ public:
 	std::vector<size_t> bounds;
 	std::vector<ValueType> entries;
 
-	inline size_t size() const { return bounds.size(); }
+	size_t size() const { return bounds.size(); }
 
-	inline typename std::vector<ValueType>::const_iterator cbegin(size_t index) const {
+	typename std::vector<ValueType>::const_iterator cbegin(size_t index) const {
 		assert(index < size());
 		return index == 0 ? entries.cbegin() : entries.cbegin() + bounds[index - 1];
 	}
 
-	inline typename std::vector<ValueType>::const_iterator cend(size_t index) const {
+	typename std::vector<ValueType>::const_iterator cend(size_t index) const {
 		assert(index < size());
 		return entries.cbegin() + bounds[index];
 	}
 
-	template <typename Iterator> inline void append(Iterator begin, Iterator end) {
+	template <typename Iterator> void append(Iterator begin, Iterator end) {
 		for (Iterator it = begin; it != end; ++it) { entries.push_back(*it); }
 		bounds.push_back(entries.size());
 	}
 
-	inline void append() { bounds.push_back(entries.size()); }
+	void append() { bounds.push_back(entries.size()); }
 
-	inline void push_back(ValueType e) {
+	void push_back(ValueType e) {
 		assert(0 < size());
 		entries.push_back(e);
 		++bounds.back();
 	}
 
-	inline void pop_back() {
+	void pop_back() {
 		assert(0 < size());
 		entries.pop_back();
 		--bounds.back();
 	}
 
-	template <typename Collection> inline void append(const Collection collection) {
+	template <typename Collection> void append(const Collection collection) {
 		append(collection.cbegin(), collection.cend());
 	}
 };
 
 template <typename Heap>
-inline void push_entry(Heap& column, index_t i, coefficient_t c, value_t diameter) {
+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));
 }
@@ -514,7 +493,7 @@ void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce,
 		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(columns_to_reduce[i]);
+		value_t diameter = get_diameter(column_to_reduce);
 
 #ifdef INDICATE_PROGRESS
 		std::cout << "\033[K"
@@ -523,11 +502,11 @@ void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce,
 #endif
 
 		index_t j = i;
-		auto column_to_add = column_to_reduce;
+		auto& column_to_add = column_to_reduce;
 
 		// start with a dummy pivot entry with coefficient -1 in order to initialize
 		// working_coboundary with the coboundary of the simplex with index column_to_reduce
-		diameter_entry_t pivot = make_diameter_entry(0, -1, modulus - 1);
+		diameter_entry_t pivot = make_diameter_entry(0, -1, -1 + modulus);
 
 #ifdef ASSEMBLE_REDUCTION_MATRIX
 		reduction_matrix.append();
@@ -552,7 +531,7 @@ void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce,
 #endif
 			{
 #ifdef ASSEMBLE_REDUCTION_MATRIX
-				auto simplex = *it;
+				const auto& simplex = *it;
 				coefficient_t simplex_coefficient = get_coefficient(simplex);
 				simplex_coefficient *= factor;
 				simplex_coefficient %= modulus;
@@ -560,7 +539,7 @@ void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce,
 #ifdef USE_COEFFICIENTS
 				const auto& simplex = reduction_coefficients[j];
 #else
-				const auto simplex = column_to_add;
+				const auto& simplex = column_to_add;
 #endif
 #endif
 				value_t simplex_diameter = get_diameter(simplex);
@@ -602,7 +581,7 @@ void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce,
 						if (might_be_apparent_pair && (simplex_diameter == coface_diameter)) {
 							if (pivot_column_index.find(coface_index) == pivot_column_index.end()) {
 								pivot = coface_entry;
-								goto found_pivot;
+								goto found_persistence_pair;
 							}
 							might_be_apparent_pair = false;
 						}
@@ -613,68 +592,66 @@ void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce,
 
 			pivot = get_pivot(working_coboundary, modulus);
 
-		found_pivot:
 			if (get_index(pivot) != -1) {
 				auto pair = pivot_column_index.find(get_index(pivot));
 
-				if (pair == pivot_column_index.end()) {
-					pivot_column_index.insert(std::make_pair(get_index(pivot), i));
+				if (pair != pivot_column_index.end()) {
+					j = pair->second;
+					column_to_add = columns_to_reduce[j];
+					continue;
+				}
+			} else {
+#ifdef PRINT_PERSISTENCE_PAIRS
+#ifdef INDICATE_PROGRESS
+				std::cout << "\033[K";
+#endif
+				std::cout << " [" << diameter << ", )" << std::endl << std::flush;
+#endif
+				break;
+			}
+
+		found_persistence_pair:
+#ifdef PRINT_PERSISTENCE_PAIRS
+			value_t death = get_diameter(pivot);
+			if (diameter != death) {
+#ifdef INDICATE_PROGRESS
+				std::cout << "\033[K";
+#endif
+				std::cout << " [" << diameter << "," << death << ")" << std::endl << std::flush;
+			}
+#endif
+
+			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_matrix (with a single diagonal 1 entry)
-					// by reduction_column (possibly with a different entry on the diagonal)
-					reduction_matrix.pop_back();
-					while (true) {
-						diameter_entry_t e = pop_pivot(reduction_column, modulus);
-						index_t index = get_index(e);
-						if (index == -1) break;
+			// replace current column of reduction_matrix (with a single diagonal 1 entry)
+			// by reduction_column (possibly with a different entry on the diagonal)
+			reduction_matrix.pop_back();
+			while (true) {
+				diameter_entry_t e = pop_pivot(reduction_column, modulus);
+				index_t index = get_index(e);
+				if (index == -1) break;
 #ifdef USE_COEFFICIENTS
-						const coefficient_t coefficient = inverse * get_coefficient(e) % modulus;
-						assert(coefficient > 0);
+				const coefficient_t coefficient = inverse * get_coefficient(e) % modulus;
+				assert(coefficient > 0);
 #else
-						const coefficient_t coefficient = 1;
+				const coefficient_t coefficient = 1;
 #endif
-						reduction_matrix.push_back(
-						    make_diameter_entry(get_diameter(e), index, coefficient));
-					}
+				reduction_matrix.push_back(
+				    make_diameter_entry(get_diameter(e), index, coefficient));
+			}
 #else
 #ifdef USE_COEFFICIENTS
-					reduction_coefficients.pop_back();
-					reduction_coefficients.push_back(
-					    make_diameter_entry(column_to_reduce, inverse));
-#endif
-#endif
-
-#ifdef PRINT_PERSISTENCE_PAIRS
-					value_t death = get_diameter(pivot);
-					if (diameter != death) {
-#ifdef INDICATE_PROGRESS
-						std::cout << "\033[K";
+			reduction_coefficients.pop_back();
+			reduction_coefficients.push_back(make_diameter_entry(column_to_reduce, inverse));
 #endif
-						std::cout << " [" << diameter << "," << death << ")" << std::endl
-						          << std::flush;
-					}
-#endif
-					break;
-				}
-				j = pair->second;
-				column_to_add = columns_to_reduce[j];
-			}
-		} while (get_index(pivot) != -1);
-
-#ifdef PRINT_PERSISTENCE_PAIRS
-		if (get_index(pivot) == -1) {
-			value_t birth = diameter;
-#ifdef INDICATE_PROGRESS
-			std::cout << "\033[K";
-#endif
-			std::cout << " [" << birth << ", )" << std::endl << std::flush;
-		}
 #endif
+			break;
+		} while (true);
 	}
 
 #ifdef INDICATE_PROGRESS
@@ -758,7 +735,6 @@ int main(int argc, char** argv) {
 	std::vector<value_t> diameters;
 
 #ifdef FILE_FORMAT_DIPHA
-
 	int64_t magic_number;
 	input_stream.read(reinterpret_cast<char*>(&magic_number), sizeof(int64_t));
 	if (magic_number != 8067171840) {
@@ -790,11 +766,9 @@ int main(int argc, char** argv) {
 
 	compressed_lower_distance_matrix_adapter dist(diameters, dist_full);
 	std::cout << "distance matrix transformed to lower triangular form" << std::endl;
-
 #endif
 
 #ifdef FILE_FORMAT_UPPER_TRIANGULAR_CSV
-
 	std::vector<value_t> distances;
 	std::string value_string;
 	while (std::getline(input_stream, value_string, ','))
@@ -807,11 +781,9 @@ int main(int argc, char** argv) {
 
 	compressed_lower_distance_matrix_adapter dist(diameters, dist_upper);
 	std::cout << "distance matrix transformed to lower triangular form" << std::endl;
-
 #endif
 
 #ifdef FILE_FORMAT_LOWER_TRIANGULAR_CSV
-
 	std::vector<value_t>& distances = diameters;
 	std::string value_string;
 	while (std::getline(input_stream, value_string, ','))
@@ -822,7 +794,6 @@ int main(int argc, char** argv) {
 	index_t n = dist.size();
 
 	std::cout << "distance matrix with " << n << " points" << std::endl;
-
 #endif
 
 	auto result = std::minmax_element(dist.distances.begin(), dist.distances.end());