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

2 files changed, 146 insertions, 38 deletions

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<typename Filtration_value>
 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<typename InputPointRange, typename Point_d >
-  Rips_complex(const InputPointRange& points, Filtration_value threshold,
-               Filtration_value distance(const Point_d& p1,const Point_d& p2)) {
+  template<typename InputPointRange, typename Distance >
+  Rips_complex(const InputPointRange& points, Filtration_value threshold, Distance distance) {
+    compute_proximity_graph<InputPointRange, Distance >(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 <typename SimplicialComplexForRips>
+  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 <typename SimplicialComplexForRips>
-  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_