Merge last trunk modifications

git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/cubical_complex_small_fix@4088 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: 422db01787f7c105092e7f12b023a7a8e1f7d316

11 files changed, 222 insertions, 48 deletions

diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h
index 30f89d65..c3085dff 100644
--- a/src/Nerve_GIC/include/gudhi/GIC.h
+++ b/src/Nerve_GIC/include/gudhi/GIC.h
@@ -230,18 +230,17 @@ class Cover_complex {
 
   /** \brief Reads and stores the input point cloud from vector stored in memory.
    *
-   * @param[in] cloud input vector representing the point cloud. Each row is a point and each coordinate is a vector.
+   * @param[in] point_cloud input vector representing the point cloud. Each row is a point and each coordinate is a vector.
    *
    */
-  template <class InputRange>
-  void set_point_cloud_from_range(InputRange const & cloud) {
-    n = cloud.size(); data_dimension = cloud[0].size(); point_cloud_name = "matrix";
+  void set_point_cloud_from_range(const std::vector<std::vector<double> > & point_cloud) {
+    n = point_cloud.size(); data_dimension = point_cloud[0].size();
+    point_cloud_name = "matrix"; cover.resize(n);
     for(int i = 0; i < n; i++){
-      point_cloud.emplace_back(cloud[i].begin(), cloud[i].begin() + data_dimension);
       boost::add_vertex(one_skeleton_OFF);
       vertices.push_back(boost::add_vertex(one_skeleton));
-      cover.emplace_back();
     }
+    this->point_cloud = point_cloud;
   }
 
   /** \brief Reads and stores the input point cloud from .(n)OFF file.
@@ -395,18 +394,12 @@ class Cover_complex {
    * @param[in] distance_matrix input vector representing the distance matrix.
    *
    */
-  template <class InputRange>
-  void set_distances_from_range(InputRange const & distance_matrix) {
-    if(point_cloud.size() == 0){
-      n = distance_matrix.size();
-      point_cloud_name = "matrix";
-      data_dimension = 0;
-      for(int i = 0; i < n; i++){
-        point_cloud.emplace_back();
-        boost::add_vertex(one_skeleton_OFF);
-        vertices.push_back(boost::add_vertex(one_skeleton));
-        cover.emplace_back();
-      }
+  void set_distances_from_range(const std::vector<std::vector<double> > & distance_matrix) {
+    n = distance_matrix.size(); data_dimension = 0; point_cloud_name = "matrix";
+    cover.resize(n); point_cloud.resize(n);
+    for(int i = 0; i < n; i++){
+      boost::add_vertex(one_skeleton_OFF);
+      vertices.push_back(boost::add_vertex(one_skeleton));
     }
     distances = distance_matrix;
   }
diff --git a/src/Persistent_cohomology/example/persistence_from_file.cpp b/src/Persistent_cohomology/example/persistence_from_file.cpp
index 53456919..0a05c193 100644
--- a/src/Persistent_cohomology/example/persistence_from_file.cpp
+++ b/src/Persistent_cohomology/example/persistence_from_file.cpp
@@ -33,7 +33,6 @@
 using namespace Gudhi;
 using namespace Gudhi::persistent_cohomology;
 
-typedef int Vertex_handle;
 typedef double Filtration_value;
 
 void program_options(int argc, char * argv[]
diff --git a/src/Rips_complex/include/gudhi/Rips_complex.h b/src/Rips_complex/include/gudhi/Rips_complex.h
index f0fe57f4..e902e52c 100644
--- a/src/Rips_complex/include/gudhi/Rips_complex.h
+++ b/src/Rips_complex/include/gudhi/Rips_complex.h
@@ -59,7 +59,7 @@ class Rips_complex {
   /**
    * \brief Type of the one skeleton graph stored inside the Rips complex structure.
    */
-  typedef typename boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS
+  typedef typename boost::adjacency_list < boost::vecS, boost::vecS, boost::directedS
   , boost::property < vertex_filtration_t, Filtration_value >
   , boost::property < edge_filtration_t, Filtration_value >> OneSkeletonGraph;
 
diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h
index 4dcc08ed..00da148f 100644
--- a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h
+++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h
@@ -55,7 +55,7 @@ template <typename Filtration_value>
 class Sparse_rips_complex {
  private:
   // TODO(MG): use a different graph where we know we can safely insert in parallel.
-  typedef typename boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS,
+  typedef typename boost::adjacency_list<boost::vecS, boost::vecS, boost::directedS,
                                          boost::property<vertex_filtration_t, Filtration_value>,
                                          boost::property<edge_filtration_t, Filtration_value>>
       Graph;
diff --git a/src/Simplex_tree/example/cech_complex_cgal_mini_sphere_3d.cpp b/src/Simplex_tree/example/cech_complex_cgal_mini_sphere_3d.cpp
index 34092ef6..6bab8adb 100644
--- a/src/Simplex_tree/example/cech_complex_cgal_mini_sphere_3d.cpp
+++ b/src/Simplex_tree/example/cech_complex_cgal_mini_sphere_3d.cpp
@@ -50,7 +50,7 @@ using Vertex_handle = Simplex_tree::Vertex_handle;
 using Simplex_handle = Simplex_tree::Simplex_handle;
 using Filtration_value = Simplex_tree::Filtration_value;
 using Siblings = Simplex_tree::Siblings;
-using Graph_t = boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS,
+using Graph_t = boost::adjacency_list<boost::vecS, boost::vecS, boost::directedS,
                                       boost::property<Gudhi::vertex_filtration_t, Filtration_value>,
                                       boost::property<Gudhi::edge_filtration_t, Filtration_value> >;
 using Edge_t = std::pair<Vertex_handle, Vertex_handle>;
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h
index 78697ea2..4b18651c 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h
@@ -1049,7 +1049,8 @@ class Simplex_tree {
    * boost::graph_traits<OneSkeletonGraph>::vertex_descriptor
    *                                    must be Vertex_handle.
    * boost::graph_traits<OneSkeletonGraph>::directed_category
-   *                                    must be undirected_tag.
+   *                                    can be directed_tag (the fastest, the least RAM use), undirected_tag or even
+   *                                    bidirected_tag.
    *
    * If an edge appears with multiplicity, the function will arbitrarily pick
    * one representative to read the filtration value.  */
@@ -1077,12 +1078,14 @@ class Simplex_tree {
                                   root_.members_.end(), *v_it,
                                   Node(&root_, boost::get(vertex_filtration_t(), skel_graph, *v_it)));
     }
-    typename boost::graph_traits<OneSkeletonGraph>::edge_iterator e_it,
-        e_it_end;
-    for (std::tie(e_it, e_it_end) = boost::edges(skel_graph); e_it != e_it_end;
-         ++e_it) {
-      auto u = source(*e_it, skel_graph);
-      auto v = target(*e_it, skel_graph);
+    std::pair<typename boost::graph_traits<OneSkeletonGraph>::edge_iterator,
+              typename boost::graph_traits<OneSkeletonGraph>::edge_iterator> boost_edges = boost::edges(skel_graph);
+    // boost_edges.first is the equivalent to boost_edges.begin()
+    // boost_edges.second is the equivalent to boost_edges.end()
+    for (; boost_edges.first != boost_edges.second; boost_edges.first++) {
+      auto edge = *(boost_edges.first);
+      auto u = source(edge, skel_graph);
+      auto v = target(edge, skel_graph);
       if (u == v) throw "Self-loops are not simplicial";
       // We cannot skip edges with the wrong orientation and expect them to
       // come a second time with the right orientation, that does not always
@@ -1098,7 +1101,7 @@ class Simplex_tree {
       }
 
       sh->second.children()->members().emplace(v,
-          Node(sh->second.children(), boost::get(edge_filtration_t(), skel_graph, *e_it)));
+          Node(sh->second.children(), boost::get(edge_filtration_t(), skel_graph, edge)));
     }
   }
 
diff --git a/src/Simplex_tree/test/simplex_tree_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_unit_test.cpp
index f63ea080..b27bbce8 100644
--- a/src/Simplex_tree/test/simplex_tree_unit_test.cpp
+++ b/src/Simplex_tree/test/simplex_tree_unit_test.cpp
@@ -864,34 +864,60 @@ BOOST_AUTO_TEST_CASE(make_filtration_non_decreasing) {
   
 }
 
-BOOST_AUTO_TEST_CASE(insert_graph) {
+
+typedef boost::mpl::list<boost::adjacency_list<boost::setS, boost::vecS, boost::directedS,
+                                               boost::property<vertex_filtration_t, double>,
+                                               boost::property<edge_filtration_t, double>>,
+                         boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS,
+                                               boost::property<vertex_filtration_t, double>,
+                                               boost::property<edge_filtration_t, double>>,
+                         boost::adjacency_list<boost::setS, boost::vecS, boost::bidirectionalS,
+                                               boost::property<vertex_filtration_t, double>,
+                                               boost::property<edge_filtration_t, double>>,
+                         boost::adjacency_list<boost::vecS, boost::vecS, boost::directedS,
+                                               boost::property<vertex_filtration_t, double>,
+                                               boost::property<edge_filtration_t, double>>,
+                         boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS,
+                                               boost::property<vertex_filtration_t, double>,
+                                               boost::property<edge_filtration_t, double>>,
+                         boost::adjacency_list<boost::vecS, boost::vecS, boost::bidirectionalS,
+                                               boost::property<vertex_filtration_t, double>,
+                                               boost::property<edge_filtration_t, double>>> list_of_graph_variants;
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(simplex_tree_insert_graph, Graph, list_of_graph_variants) {
   std::cout << "********************************************************************" << std::endl;
   std::cout << "INSERT GRAPH" << std::endl;
-  typedef typename boost::adjacency_list<boost::vecS, boost::vecS,
-          boost::undirectedS,
-          boost::property<vertex_filtration_t, double>,
-          boost::property<edge_filtration_t, double>> Graph;
+
   Graph g(3);
   // filtration value 0 everywhere
   put(Gudhi::vertex_filtration_t(), g, 0, 0);
   put(Gudhi::vertex_filtration_t(), g, 1, 0);
   put(Gudhi::vertex_filtration_t(), g, 2, 0);
   // vertices don't always occur in sorted order
-  add_edge(0, 1, 0, g);
-  add_edge(2, 1, 0, g);
-  add_edge(2, 0, 0, g);
+  add_edge(0, 1, 1.1, g);
+  add_edge(2, 0, 2.2, g);
+  add_edge(2, 1, 3.3, g);
 
-  typedef Simplex_tree<> typeST;
-  typeST st1;
+  Simplex_tree<> st1;
   st1.insert_graph(g);
   BOOST_CHECK(st1.num_simplices() == 6);
 
   // edges can have multiplicity in the graph unless we replace the first vecS with (hash_)setS
-  add_edge(1, 0, 0, g);
-  add_edge(1, 2, 0, g);
-  add_edge(0, 2, 0, g);
-  add_edge(0, 2, 0, g);
-  typeST st2;
+  add_edge(1, 0, 1.1, g);
+  add_edge(1, 2, 3.3, g);
+  add_edge(0, 2, 2.2, g);
+  add_edge(0, 1, 1.1, g);
+  add_edge(2, 1, 3.3, g);
+  add_edge(2, 0, 2.2, g);
+  Simplex_tree<> st2;
   st2.insert_graph(g);
   BOOST_CHECK(st2.num_simplices() == 6);
+
+  std::cout << "st1 is" << std::endl;
+  std::cout << st1 << std::endl;
+
+  std::cout << "st2 is" << std::endl;
+  std::cout << st2 << std::endl;
+
+  BOOST_CHECK(st1 == st2);
 }
diff --git a/src/common/benchmark/CMakeLists.txt b/src/common/benchmark/CMakeLists.txt
new file mode 100644
index 00000000..a3787d6e
--- /dev/null
+++ b/src/common/benchmark/CMakeLists.txt
@@ -0,0 +1,3 @@
+project(common_benchmark)
+
+add_executable(Graph_simplicial_complex_benchmark Graph_simplicial_complex_benchmark.cpp)
diff --git a/src/common/benchmark/Graph_simplicial_complex_benchmark.cpp b/src/common/benchmark/Graph_simplicial_complex_benchmark.cpp
new file mode 100644
index 00000000..825d6cb5
--- /dev/null
+++ b/src/common/benchmark/Graph_simplicial_complex_benchmark.cpp
@@ -0,0 +1,150 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2018 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
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <gudhi/graph_simplicial_complex.h>
+#include <gudhi/distance_functions.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Clock.h>
+#include <gudhi/Points_off_io.h>
+
+#include <boost/graph/adjacency_list.hpp>
+
+#include <iostream>
+#include <string>
+#include <vector>
+#include <limits>  // for numeric limits
+#include <fstream>
+#include <cassert>
+
+
+std::ofstream results_csv("results.csv");
+
+template< typename Adjacency_list, typename ForwardPointRange, typename Distance >
+Adjacency_list proximity_graph_computation(const ForwardPointRange& points, double threshold, Distance distance) {
+  std::vector<std::pair< int, int >> edges;
+  std::vector< double > edges_fil;
+  std::map< int, double > vertices;
+
+  int idx_u, idx_v;
+  double fil;
+  idx_u = 0;
+  for (auto it_u = points.begin(); it_u != points.end(); ++it_u) {
+    idx_v = idx_u + 1;
+    for (auto it_v = it_u + 1; it_v != points.end(); ++it_v, ++idx_v) {
+      fil = distance(*it_u, *it_v);
+      if (fil <= threshold) {
+        edges.emplace_back(idx_u, idx_v);
+        edges_fil.push_back(fil);
+      }
+    }
+    ++idx_u;
+  }
+
+  // Points are labeled from 0 to idx_u-1
+  Adjacency_list skel_graph(edges.begin(), edges.end(), edges_fil.begin(), idx_u);
+
+  auto vertex_prop = boost::get(Gudhi::vertex_filtration_t(), skel_graph);
+
+  typename boost::graph_traits<Adjacency_list>::vertex_iterator vi, vi_end;
+  for (std::tie(vi, vi_end) = boost::vertices(skel_graph);
+       vi != vi_end; ++vi) {
+    boost::put(vertex_prop, *vi, 0.);
+  }
+
+  return skel_graph;
+}
+
+template <typename Adjacency_list>
+void benchmark_proximity_graph(const std::string& msg, const std::string& off_file_name) {
+  Gudhi::Points_off_reader<std::vector<double>> off_reader(off_file_name);
+  assert(off_reader.is_valid());
+
+  std::cout << "+ " << msg << std::endl;
+
+  results_csv << "\"nb_points\";"
+              << "\"nb_simplices\";"
+              << "\"compute proximity graph(sec.)\";"
+              << "\"complex_creation_time(sec.)\";"
+              << "\"" << msg << "\";" << std::endl;
+
+  Gudhi::Clock pg_compute_proximity_graph("    benchmark_proximity_graph - compute proximity graph");
+  pg_compute_proximity_graph.begin();
+  // benchmark begin
+  Adjacency_list proximity_graph = proximity_graph_computation<Adjacency_list>(off_reader.get_point_cloud(),
+                                                                               std::numeric_limits<double>::infinity(),
+                                                                               Gudhi::Euclidean_distance());
+  // benchmark end
+  pg_compute_proximity_graph.end();
+  std::cout << pg_compute_proximity_graph;
+
+  Gudhi::Simplex_tree<> complex;
+  Gudhi::Clock st_create_clock("    benchmark_proximity_graph - complex creation");
+  st_create_clock.begin();
+  // benchmark begin
+  complex.insert_graph(proximity_graph);
+  // benchmark end
+  st_create_clock.end();
+  std::cout << st_create_clock;
+
+  results_csv << off_reader.get_point_cloud().size() << ";" << complex.num_simplices() << ";"
+              << pg_compute_proximity_graph.num_seconds() << ";"
+              << st_create_clock.num_seconds() << ";" << std::endl;
+
+  std::cout << "    benchmark_proximity_graph - nb simplices = " << complex.num_simplices() << std::endl;
+}
+
+int main(int argc, char * const argv[]) {
+  std::string off_file_name(argv[1]);
+
+  // The fastest, the less memory used
+  using vecSdirectedS = boost::adjacency_list<boost::vecS, boost::vecS, boost::directedS,
+                                              boost::property<Gudhi::vertex_filtration_t, double>,
+                                              boost::property<Gudhi::edge_filtration_t, double>>;
+  benchmark_proximity_graph<vecSdirectedS>("vecSdirectedS", off_file_name);
+
+  using vecSundirectedS = boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS,
+                                                boost::property<Gudhi::vertex_filtration_t, double>,
+                                                boost::property<Gudhi::edge_filtration_t, double>>;
+  benchmark_proximity_graph<vecSundirectedS>("vecSundirectedS", off_file_name);
+
+  using vecSbidirectionalS = boost::adjacency_list<boost::vecS, boost::vecS, boost::bidirectionalS,
+                                                   boost::property<Gudhi::vertex_filtration_t, double>,
+                                                   boost::property<Gudhi::edge_filtration_t, double>>;
+  benchmark_proximity_graph<vecSbidirectionalS>("vecSbidirectionalS", off_file_name);
+
+  using setSdirectedS = boost::adjacency_list<boost::setS, boost::vecS, boost::directedS,
+                                              boost::property<Gudhi::vertex_filtration_t, double>,
+                                              boost::property<Gudhi::edge_filtration_t, double>>;
+  benchmark_proximity_graph<setSdirectedS>("setSdirectedS", off_file_name);
+
+  using setSundirectedS = boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS,
+                                                boost::property<Gudhi::vertex_filtration_t, double>,
+                                                boost::property<Gudhi::edge_filtration_t, double>>;
+  benchmark_proximity_graph<setSundirectedS>("setSundirectedS", off_file_name);
+
+  using setSbidirectionalS = boost::adjacency_list<boost::setS, boost::vecS, boost::bidirectionalS,
+                                                   boost::property<Gudhi::vertex_filtration_t, double>,
+                                                   boost::property<Gudhi::edge_filtration_t, double>>;
+  benchmark_proximity_graph<setSbidirectionalS>("setSbidirectionalS", off_file_name);
+
+  return 0;
+}
diff --git a/src/common/include/gudhi/graph_simplicial_complex.h b/src/common/include/gudhi/graph_simplicial_complex.h
index 49fe56cc..0d81ca71 100644
--- a/src/common/include/gudhi/graph_simplicial_complex.h
+++ b/src/common/include/gudhi/graph_simplicial_complex.h
@@ -49,7 +49,7 @@ struct vertex_filtration_t {
  *
  */
 template <typename SimplicialComplexForProximityGraph>
-using Proximity_graph = typename boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS
+using Proximity_graph = typename boost::adjacency_list < boost::vecS, boost::vecS, boost::directedS
 , boost::property < vertex_filtration_t, typename SimplicialComplexForProximityGraph::Filtration_value >
 , boost::property < edge_filtration_t, typename SimplicialComplexForProximityGraph::Filtration_value >>;
 
diff --git a/src/common/include/gudhi/reader_utils.h b/src/common/include/gudhi/reader_utils.h
index 26eeb76d..0ee7649d 100644
--- a/src/common/include/gudhi/reader_utils.h
+++ b/src/common/include/gudhi/reader_utils.h
@@ -172,10 +172,10 @@ bool read_simplex(std::istream& in_, std::vector<Vertex_handle>& simplex, Filtra
   if (!(in_ >> dim)) return false;
   Vertex_handle v;
   for (int i = 0; i < dim + 1; ++i) {
-    in_ >> v;
+    if (!(in_ >> v)) return false;
     simplex.push_back(v);
   }
-  in_ >> fil;
+  if (!(in_ >> fil)) return false;
   in_.ignore((std::numeric_limits<std::streamsize>::max)(), '\n');  // ignore until the carriage return
   return true;
 }