point cloud support

2 files changed, 131 insertions, 143 deletions

diff --git a/ripser.cpp b/ripser.cpp
index 78035ba..7d4c28d 100644
--- a/ripser.cpp
+++ b/ripser.cpp
@@ -15,7 +15,6 @@
     You should have received a copy of the GNU General Public License
     along with this program.  If not, see <http://www.gnu.org/licenses/>. */
 
-
 //#define ASSEMBLE_REDUCTION_MATRIX
 //#define USE_COEFFICIENTS
 
@@ -30,13 +29,13 @@
 //#define USE_GOOGLE_HASHMAP
 
 #include <algorithm>
-#include <iostream>
-#include <fstream>
-#include <sstream>
 #include <cassert>
-#include <queue>
 #include <cmath>
+#include <fstream>
+#include <iostream>
 #include <numeric>
+#include <queue>
+#include <sstream>
 #include <unordered_map>
 
 #ifdef USE_GOOGLE_HASHMAP
@@ -301,26 +300,25 @@ public:
 		init_rows();
 	}
 
-    value_t operator()(const index_t i, const index_t j) const;
+	value_t operator()(const index_t i, const index_t j) const;
 
 	size_t size() const { return rows.size(); }
 };
 
 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)) { }
-    
-    value_t operator()(const index_t i, const index_t j) const {
-        value_t result;
-        result = std::inner_product((points[i]).begin(), (points[i]).end(), (points[j]).begin(), result,
-                                    std::plus<value_t>(), [] (value_t u, value_t v) { return (u-v)*(u-v); });
-        return std::sqrt(result);
-    }
-    
-    size_t size() const { return points.size(); }
+	std::vector<std::vector<value_t>> 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); }));
+	}
+
+	size_t size() const { return points.size(); }
 };
 
 template <> void compressed_distance_matrix<LOWER_TRIANGULAR>::init_rows() {
@@ -341,7 +339,7 @@ template <> void compressed_distance_matrix<UPPER_TRIANGULAR>::init_rows() {
 
 template <>
 value_t compressed_distance_matrix<UPPER_TRIANGULAR>::operator()(const index_t i,
-                                                                         const index_t j) const {
+                                                                 const index_t j) const {
 	if (i < j)
 		return rows[i][j];
 	else if (i > j)
@@ -352,7 +350,7 @@ value_t compressed_distance_matrix<UPPER_TRIANGULAR>::operator()(const index_t i
 
 template <>
 value_t compressed_distance_matrix<LOWER_TRIANGULAR>::operator()(const index_t i,
-                                                                         const index_t j) const {
+                                                                 const index_t j) const {
 	if (i > j)
 		return rows[i][j];
 	else if (i < j)
@@ -361,10 +359,8 @@ value_t compressed_distance_matrix<LOWER_TRIANGULAR>::operator()(const index_t i
 		return 0;
 }
 
-typedef compressed_distance_matrix<LOWER_TRIANGULAR>
-    compressed_lower_distance_matrix;
-typedef compressed_distance_matrix<UPPER_TRIANGULAR>
-    compressed_upper_distance_matrix;
+typedef compressed_distance_matrix<LOWER_TRIANGULAR> compressed_lower_distance_matrix;
+typedef compressed_distance_matrix<UPPER_TRIANGULAR> compressed_upper_distance_matrix;
 
 template <typename Heap> diameter_entry_t pop_pivot(Heap& column, coefficient_t modulus) {
 	if (column.empty())
@@ -460,31 +456,31 @@ void assemble_columns_to_reduce(std::vector<diameter_index_t>& columns_to_reduce
                                 hash_map<index_t, index_t>& pivot_column_index,
                                 const Comparator& comp, index_t dim, index_t n, value_t threshold,
                                 const binomial_coeff_table& binomial_coeff) {
-    index_t num_simplices = binomial_coeff(n, dim + 2);
-    
-    columns_to_reduce.clear();
-    
+	index_t num_simplices = binomial_coeff(n, dim + 2);
+
+	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 = comp.diameter(index);
-            if (diameter <= threshold) columns_to_reduce.push_back(std::make_pair(diameter, index));
-        }
-    }
-    
+	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 = comp.diameter(index);
+			if (diameter <= threshold) columns_to_reduce.push_back(std::make_pair(diameter, index));
+		}
+	}
+
 #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
 }
 
@@ -516,11 +512,13 @@ void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce,
 
 #ifdef ASSEMBLE_REDUCTION_MATRIX
 		std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>,
-		                    smaller_index<diameter_entry_t>> reduction_column;
+		                    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;
+		                    greater_diameter_or_smaller_index<diameter_entry_t>>
+		    working_coboundary;
 
 		value_t diameter = get_diameter(column_to_reduce);
 
@@ -762,111 +760,100 @@ int main(int argc, char** argv) {
 		}
 	}
 
-	std::ifstream input_stream(filename, std::ios_base::binary | std::ios_base::in);
+	std::ifstream input_stream(filename);
 	if (input_stream.fail()) {
 		std::cerr << "couldn't open file " << filename << std::endl;
 		exit(-1);
 	}
-    
+
 #ifdef FILE_FORMAT_POINT_CLOUD
-    
-    std::vector<std::vector<value_t>> points;
-
-    std::string line;
-    value_t value;
-    while (std::getline(input_stream, line)) {
-        std::vector<value_t> point;
-        std::istringstream s(line);
-        while (s >> value) point.push_back(value);
-        points.push_back(point);
-    }
-    
-//    euclidean_distance_matrix dist(points);
-//    index_t n = dist.size();
-
-    euclidean_distance_matrix eucl_dist(points);
-    
-    index_t n = eucl_dist.size();
-
-    std::vector<value_t> distances;
-    
-    for (int i = 0; i < n; ++i)
-        for (int j = 0; j < i; ++j)
-            if (i > j) distances.push_back(eucl_dist(i,j));
-    
-    compressed_lower_distance_matrix dist(distances);
-    
-    std::cout << "point cloud with " << n << " points" << std::endl;
-    
-#endif
-    
-    
+
+	std::vector<std::vector<value_t>> points;
+
+	std::string line;
+	value_t value;
+	while (std::getline(input_stream, line)) {
+		std::vector<value_t> point;
+		std::istringstream s(line);
+		while (s >> value) point.push_back(value);
+		if (!point.empty()) points.push_back(point);
+	}
+
+	euclidean_distance_matrix dist(points);
+	index_t n = dist.size();
+
+	std::cout << "point cloud with " << n << " points" << std::endl;
+
+#endif
+
 #ifdef FILE_FORMAT_LOWER_TRIANGULAR_CSV
-    
-    std::vector<value_t> distances;
-    value_t value;
-    while (input_stream >> value) {
-        distances.push_back(value);
-        input_stream.ignore();
-    }
-    
-    compressed_lower_distance_matrix dist(distances);
-    
-    index_t n = dist.size();
-    
-    std::cout << "distance matrix with " << n << " points" << std::endl;
-    
-#endif
-    
+
+	std::vector<value_t> distances;
+	value_t value;
+	while (input_stream >> value) {
+		distances.push_back(value);
+		input_stream.ignore();
+	}
+
+	compressed_lower_distance_matrix dist(distances);
+
+	index_t n = dist.size();
+
+	std::cout << "distance matrix with " << n << " points" << std::endl;
+
+#endif
+
 #ifdef FILE_FORMAT_UPPER_TRIANGULAR_CSV
-    
-    std::vector<value_t> distances;
-    value_t value;
-    while (input_stream >> value) {
-        distances.push_back(value);
-        input_stream.ignore();
-    }
-    
-    compressed_upper_distance_matrix dist(distances);
-    
-    index_t n = dist.size();
-    
-    std::cout << "distance matrix with " << n << " points" << std::endl;
-    
-#endif
-
-    
+
+	std::vector<value_t> distances;
+	value_t value;
+	while (input_stream >> value) {
+		distances.push_back(value);
+		input_stream.ignore();
+	}
+
+	compressed_upper_distance_matrix dist(distances);
+
+	index_t n = dist.size();
+
+	std::cout << "distance matrix with " << n << " points" << std::endl;
+
+#endif
+
 #ifdef FILE_FORMAT_DIPHA
-    
-    if (read<int64_t>(input_stream) != 8067171840) {
-        std::cerr << filename << " is not a Dipha file (magic number: 8067171840)" << std::endl;
-        exit(-1);
-    }
-    
-    if (read<int64_t>(input_stream) != 7) {
-        std::cerr << filename << " is not a Dipha distance matrix (file type: 7)" << std::endl;
-        exit(-1);
-    }
-    
-    index_t n = read<int64_t>(input_stream);
-    
-    std::vector<value_t> distances;
-    
-    for (int i = 0; i < n; ++i)
-        for (int j = 0; j < n; ++j)
-            if (i > j) distances.push_back(read<double>(input_stream)); else read<double>(input_stream);
-    
-    compressed_lower_distance_matrix dist(distances);
-    
-    std::cout << "distance matrix with " << n << " points" << std::endl;
-    
-#endif
-    
+
+	if (read<int64_t>(input_stream) != 8067171840) {
+		std::cerr << filename << " is not a Dipha file (magic number: 8067171840)" << std::endl;
+		exit(-1);
+	}
+
+	if (read<int64_t>(input_stream) != 7) {
+		std::cerr << filename << " is not a Dipha distance matrix (file type: 7)" << std::endl;
+		exit(-1);
+	}
+
+	index_t n = read<int64_t>(input_stream);
+
+	std::vector<value_t> distances;
+
+	for (int i = 0; i < n; ++i)
+		for (int j = 0; j < n; ++j)
+			if (i > j)
+				distances.push_back(read<double>(input_stream));
+			else
+				read<double>(input_stream);
+
+	compressed_lower_distance_matrix dist(distances);
+
+	std::cout << "distance matrix with " << n << " points" << std::endl;
+
+#endif
+
 #ifndef FILE_FORMAT_POINT_CLOUD
 	auto result = std::minmax_element(dist.distances.begin(), dist.distances.end());
-    std::cout << "value range: [" << *result.first << "," << *result.second << "]" << std::endl;
+	std::cout << "value range: [" << *result.first << "," << *result.second << "]" << std::endl;
 #endif
-    
+
 	dim_max = std::min(dim_max, n - 2);
 
 	binomial_coeff_table binomial_coeff(n, dim_max + 2);
diff --git a/ripser.xcodeproj/project.pbxproj b/ripser.xcodeproj/project.pbxproj
index a14431e..52c2d4e 100644
--- a/ripser.xcodeproj/project.pbxproj
+++ b/ripser.xcodeproj/project.pbxproj
@@ -209,7 +209,8 @@
 			buildSettings = {
 				GCC_PREPROCESSOR_DEFINITIONS = (
 					"$(inherited)",
-					FILE_FORMAT_LOWER_TRIANGULAR_CSV,
+					_FILE_FORMAT_LOWER_TRIANGULAR_CSV,
+					FILE_FORMAT_POINT_CLOUD,
 				);
 				PRODUCT_NAME = "$(TARGET_NAME)";
 			};