From c6ab4a62da2572b51032995266c109df7fe76dfd Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Fri, 18 Nov 2016 17:21:46 +0000 Subject: Add distance matrix read from csv files Distance template instead of Point_d for distance matrices git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/distance_matrix_in_rips_module@1762 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 695bead714228c57e1be1650048125aa0b841683 --- src/Rips_complex/include/gudhi/Rips_complex.h | 89 +++++++++++++++---------- src/common/include/gudhi/reader_utils.h | 95 ++++++++++++++++++++++++++- 2 files changed, 146 insertions(+), 38 deletions(-) (limited to 'src') diff --git a/src/Rips_complex/include/gudhi/Rips_complex.h b/src/Rips_complex/include/gudhi/Rips_complex.h index da755b7c..6f947f41 100644 --- a/src/Rips_complex/include/gudhi/Rips_complex.h +++ b/src/Rips_complex/include/gudhi/Rips_complex.h @@ -2,7 +2,7 @@ * (Geometric Understanding in Higher Dimensions) is a generic C++ * library for computational topology. * - * Author(s): Clément Maria, Vincent Rouvreau + * Author(s): Clément Maria, Pawel Dlotko, Vincent Rouvreau * * Copyright (C) 2016 INRIA * @@ -58,12 +58,13 @@ template class Rips_complex { private: typedef typename boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS - , boost::property < vertex_filtration_t, Filtration_value > - , boost::property < edge_filtration_t, Filtration_value >> Graph_t; - + , boost::property < vertex_filtration_t, Filtration_value > + , boost::property < edge_filtration_t, Filtration_value >> Graph_t; + typedef int Vertex_handle; - + public: + /** \brief Rips_complex constructor from a list of points. * * @param[in] points Range of points. @@ -72,9 +73,51 @@ class Rips_complex { * * The type InputPointRange must be a range for which std::begin and std::end return input iterators on a point. */ - template - Rips_complex(const InputPointRange& points, Filtration_value threshold, - Filtration_value distance(const Point_d& p1,const Point_d& p2)) { + template + Rips_complex(const InputPointRange& points, Filtration_value threshold, Distance distance) { + compute_proximity_graph(points, threshold, distance); + } + + /** \brief Initializes the simplicial complex from the 1-skeleton graph and expands it until a given maximal + * dimension. + * + * \tparam SimplicialComplexForRips must meet `SimplicialComplexForRips` concept. + * + * @param[in] complex SimplicialComplexForRips to be created. + * @param[in] dim_max graph expansion for rips until this given maximal dimension. + * + * @return true if creation succeeds, false otherwise. + * + */ + template + bool create_complex(SimplicialComplexForRips& complex, int dim_max) { + if (complex.num_vertices() > 0) { + std::cerr << "Rips_complex create_complex - complex is not empty\n"; + return false; // ----- >> + } + + // insert the proximity graph in the simplicial complex + complex.insert_graph(rips_skeleton_graph_); + // expand the graph until dimension dim_max + complex.expansion(dim_max); + + // -------------------------------------------------------------------------------------------- + return true; + } + + public: + /** \brief Output the proximity graph of the points. + * + * If points contains n elements, the proximity graph is the graph + * with n vertices, and an edge [u,v] iff the distance function between + * points u and v is smaller than threshold. + * + * The type PointCloud furnishes .begin() and .end() methods, that return + * iterators with value_type Point. + */ + template< typename InputPointRange, typename Distance > + void compute_proximity_graph(const InputPointRange& points, Filtration_value threshold, + Distance distance) { std::vector< std::pair< Vertex_handle, Vertex_handle > > edges; std::vector< Filtration_value > edges_fil; std::map< Vertex_handle, Filtration_value > vertices; @@ -102,7 +145,7 @@ class Rips_complex { // -------------------------------------------------------------------------------------------- // Creates the proximity graph from edges and sets the property with the filtration value. // Number of points is labeled from 0 to idx_u-1 - rips_skeleton_graph_ = Graph_t(edges.begin() , edges.end() , edges_fil.begin() , idx_u); + rips_skeleton_graph_ = Graph_t(edges.begin(), edges.end(), edges_fil.begin(), idx_u); auto vertex_prop = boost::get(vertex_filtration_t(), rips_skeleton_graph_); @@ -112,37 +155,11 @@ class Rips_complex { vi != vi_end; ++vi) { boost::put(vertex_prop, *vi, 0.); } - } - /** \brief Initializes the simplicial complex from the 1-skeleton graph and expands it until a given maximal - * dimension. - * - * \tparam SimplicialComplexForRips must meet `SimplicialComplexForRips` concept. - * - * @param[in] complex SimplicialComplexForRips to be created. - * @param[in] dim_max graph expansion for rips until this given maximal dimension. - * - * @return true if creation succeeds, false otherwise. - * - */ - template - bool create_complex(SimplicialComplexForRips& complex, int dim_max) { - if (complex.num_vertices() > 0) { - std::cerr << "Rips_complex create_complex - complex is not empty\n"; - return false; // ----- >> - } - - // insert the proximity graph in the simplicial complex - complex.insert_graph(rips_skeleton_graph_); - // expand the graph until dimension dim_max - complex.expansion(dim_max); - - // -------------------------------------------------------------------------------------------- - return true; - } private: Graph_t rips_skeleton_graph_; + }; } // namespace rips_complex diff --git a/src/common/include/gudhi/reader_utils.h b/src/common/include/gudhi/reader_utils.h index 4154acc9..5897030f 100644 --- a/src/common/include/gudhi/reader_utils.h +++ b/src/common/include/gudhi/reader_utils.h @@ -2,9 +2,9 @@ * (Geometric Understanding in Higher Dimensions) is a generic C++ * library for computational topology. * - * Author(s): Clément Maria + * Author(s): Clement Maria, Pawel Dlotko * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * Copyright (C) 2014 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -196,4 +196,95 @@ bool read_hasse_simplex(std::istream & in_, std::vector< Simplex_key > & boundar return true; } +/** + * \brief Read a lower triangular distance matrix from a csv file. We assume that the .csv store the whole + * (square) matrix. + * + * \author Pawel Dlotko + * + * Square matrix file format: + * 0;D12;...;D1j + * D21;0;...;D2j + * ... + * Dj1;Dj2;...;0 + * + * lower matrix file format: + * 0 + * D21; + * D31;D32; + * ... + * Dj1;Dj2;...;Dj(j-1); + * + **/ +template< typename Filtration_value > +std::vector< std::vector< Filtration_value > > read_lower_triangular_matrix_from_csv_file(std::string filename, + const char separator=';') { + bool dbg = true; + if (dbg) { + std::cerr << "Using procedure read_lower_triangular_matrix_from_csv_file \n"; + } + std::vector< std::vector< Filtration_value > > result; + std::ifstream in; + in.open(filename.c_str()); + if (!in.is_open()) { + return result; + } + + std::string line; + + // the first line is emtpy, so we ignore it: + std::getline(in, line); + std::vector< Filtration_value > values_in_this_line; + result.push_back(values_in_this_line); + + int number_of_line = 0; + + // first, read the file line by line to a string: + while (std::getline(in, line)) { + // if line is empty, break + if (line.size() == 0) + break; + + // if the last element of a string is comma: + if (line[ line.size() - 1 ] == separator) { + // then shrink the string by one + line.pop_back(); + } + + // replace all commas with spaces + std::replace(line.begin(), line.end(), separator, ' '); + + // put the new line to a stream + std::istringstream iss(line); + // and now read the doubles. + + int number_of_entry = 0; + std::vector< Filtration_value > values_in_this_line; + while (iss.good()) { + double entry; + iss >> entry; + if (number_of_entry <= number_of_line) { + values_in_this_line.push_back(entry); + } + ++number_of_entry; + } + if (!values_in_this_line.empty())result.push_back(values_in_this_line); + ++number_of_line; + + } + in.close(); + + if (dbg) { + std::cerr << "Here is the matrix we read : \n"; + for (size_t i = 0; i != result.size(); ++i) { + for (size_t j = 0; j != result[i].size(); ++j) { + std::cerr << result[i][j] << " "; + } + std::cerr << std::endl; + } + } + + return result; +} // read_lower_triangular_matrix_from_csv_file + #endif // READER_UTILS_H_ -- cgit v1.2.3