Merge last Nerve_GIC development from eponym branch

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

Former-commit-id: c57e5f6f364b6162f607edd07f3a31577e640fcb

2 files changed, 157 insertions, 100 deletions

diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h
index 40ff7a4a..f5b67be6 100644
--- a/src/Nerve_GIC/include/gudhi/GIC.h
+++ b/src/Nerve_GIC/include/gudhi/GIC.h
@@ -23,6 +23,11 @@
 #ifndef GIC_H_
 #define GIC_H_
 
+#ifdef GUDHI_USE_TBB
+#include <tbb/parallel_for.h>
+#include <tbb/mutex.h>
+#endif
+
 #include <gudhi/Debug_utils.h>
 #include <gudhi/graph_simplicial_complex.h>
 #include <gudhi/reader_utils.h>
@@ -99,9 +104,8 @@ class Cover_complex {
   int data_dimension;                           // dimension of input data.
   int n;                                        // number of points.
 
-  std::map<int, double> func;  // function used to compute the output simplicial complex.
-  std::map<int, double>
-      func_color;                 // function used to compute the colors of the nodes of the output simplicial complex.
+  std::vector<double> func;       // function used to compute the output simplicial complex.
+  std::vector<double> func_color; // function used to compute the colors of the nodes of the output simplicial complex.
   bool functional_cover = false;  // whether we use a cover with preimages of a function or not.
 
   Graph one_skeleton_OFF;          // one-skeleton given by the input OFF file (if it exists).
@@ -114,8 +118,8 @@ class Cover_complex {
   Persistence_diagram PD;
   std::vector<double> distribution;
 
-  std::map<int, std::vector<int> >
-      cover;  // function associating to each data point its vectors of cover elements to which it belongs.
+  std::vector<std::vector<int> >
+      cover;  // function associating to each data point the vector of cover elements to which it belongs.
   std::map<int, std::vector<int> >
       cover_back;  // inverse of cover, in order to get the data points associated to a specific cover element.
   std::map<int, double> cover_std;  // standard function (induced by func) used to compute the extended persistence
@@ -140,8 +144,8 @@ class Cover_complex {
 
   // Point comparator
   struct Less {
-    Less(std::map<int, double> func) { Fct = func; }
-    std::map<int, double> Fct;
+    Less(std::vector<double> func) { Fct = func; }
+    std::vector<double> Fct;
     bool operator()(int a, int b) {
       if (Fct[a] == Fct[b])
         return a < b;
@@ -276,6 +280,7 @@ class Cover_complex {
         point_cloud.emplace_back(point.begin(), point.begin() + data_dimension);
         boost::add_vertex(one_skeleton_OFF);
         vertices.push_back(boost::add_vertex(one_skeleton));
+        std::vector<int> dummy; dummy.clear(); cover.push_back(dummy);
         i++;
       }
     }
@@ -430,17 +435,29 @@ class Cover_complex {
 
     if (distances.size() == 0) compute_pairwise_distances(distance);
 
-    // #pragma omp parallel for
-    for (int i = 0; i < N; i++) {
-      SampleWithoutReplacement(n, m, samples);
-      double hausdorff_dist = 0;
-      for (int j = 0; j < n; j++) {
-        double mj = distances[j][samples[0]];
-        for (int k = 1; k < m; k++) mj = std::min(mj, distances[j][samples[k]]);
-        hausdorff_dist = std::max(hausdorff_dist, mj);
+    #ifdef GUDHI_USE_TBB
+      tbb::parallel_for(0, N, [&](int i){
+        SampleWithoutReplacement(n, m, samples);
+        double hausdorff_dist = 0;
+        for (int j = 0; j < n; j++) {
+          double mj = distances[j][samples[0]];
+          for (int k = 1; k < m; k++) mj = std::min(mj, distances[j][samples[k]]);
+          hausdorff_dist = std::max(hausdorff_dist, mj);
+        }
+        delta += hausdorff_dist / N;
+      });
+    #else
+      for (int i = 0; i < N; i++) {
+        SampleWithoutReplacement(n, m, samples);
+        double hausdorff_dist = 0;
+        for (int j = 0; j < n; j++) {
+          double mj = distances[j][samples[0]];
+          for (int k = 1; k < m; k++) mj = std::min(mj, distances[j][samples[k]]);
+          hausdorff_dist = std::max(hausdorff_dist, mj);
+        }
+        delta += hausdorff_dist / N;
       }
-      delta += hausdorff_dist / N;
-    }
+    #endif
 
     if (verbose) std::cout << "delta = " << delta << std::endl;
     set_graph_from_rips(delta, distance);
@@ -465,7 +482,7 @@ class Cover_complex {
     while (std::getline(input, line)) {
       std::stringstream stream(line);
       stream >> f;
-      func.emplace(i, f);
+      func.push_back(f);
       i++;
     }
     functional_cover = true;
@@ -479,7 +496,7 @@ class Cover_complex {
            *
            */
   void set_function_from_coordinate(int k) {
-    for (int i = 0; i < n; i++) func.emplace(i, point_cloud[i][k]);
+    for (int i = 0; i < n; i++) func.push_back(point_cloud[i][k]);
     functional_cover = true;
     cover_name = "coordinate " + std::to_string(k);
   }
@@ -492,7 +509,7 @@ class Cover_complex {
            */
   template <class InputRange>
   void set_function_from_range(InputRange const& function) {
-    for (int i = 0; i < n; i++) func.emplace(i, function[i]);
+    for (int i = 0; i < n; i++) func.push_back(function[i]);
     functional_cover = true;
   }
 
@@ -710,37 +727,69 @@ class Cover_complex {
       funcstd[i] = 0.5 * (u + v);
     }
 
-    if (verbose) std::cout << "Computing connected components..." << std::endl;
-    // #pragma omp parallel for
-    for (int i = 0; i < res; i++) {
-      // Compute connected components
-      Graph G = one_skeleton.create_subgraph();
-      int num = preimages[i].size();
-      std::vector<int> component(num);
-      for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G);
-      boost::connected_components(G, &component[0]);
-      int max = 0;
-
-      // For each point in preimage
-      for (int j = 0; j < num; j++) {
-        // Update number of components in preimage
-        if (component[j] > max) max = component[j];
-
-        // Identify component with Cantor polynomial N^2 -> N
-        int identifier = (std::pow(i + component[j], 2) + 3 * i + component[j]) / 2;
-
-        // Update covers
-        cover[preimages[i][j]].push_back(identifier);
-        cover_back[identifier].push_back(preimages[i][j]);
-        cover_fct[identifier] = i;
-        cover_std[identifier] = funcstd[i];
-        cover_color[identifier].second += func_color[preimages[i][j]];
-        cover_color[identifier].first += 1;
-      }
+    #ifdef GUDHI_USE_TBB
+      if (verbose) std::cout << "Computing connected components (parallelized)..." << std::endl;
+      tbb::parallel_for(0, res, [&](int i){
+        // Compute connected components
+        Graph G = one_skeleton.create_subgraph();
+        int num = preimages[i].size();
+        std::vector<int> component(num);
+        for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G);
+        boost::connected_components(G, &component[0]);
+        int max = 0;
+
+        // For each point in preimage
+        for (int j = 0; j < num; j++) {
+          // Update number of components in preimage
+          if (component[j] > max) max = component[j];
+
+          // Identify component with Cantor polynomial N^2 -> N
+          int identifier = ((i + component[j])*(i + component[j]) + 3 * i + component[j]) / 2;
+
+          // Update covers
+          cover[preimages[i][j]].push_back(identifier);
+          cover_back[identifier].push_back(preimages[i][j]);
+          cover_fct[identifier] = i;
+          cover_std[identifier] = funcstd[i];
+          cover_color[identifier].second += func_color[preimages[i][j]];
+          cover_color[identifier].first += 1;
+        }
 
-      // Maximal dimension is total number of connected components
-      id += max + 1;
-    }
+        // Maximal dimension is total number of connected components
+        id += max + 1;
+      });
+    #else
+      if (verbose) std::cout << "Computing connected components..." << std::endl;
+      for (int i = 0; i < res; i++) {
+        // Compute connected components
+        Graph G = one_skeleton.create_subgraph();
+        int num = preimages[i].size();
+        std::vector<int> component(num);
+        for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G);
+        boost::connected_components(G, &component[0]);
+        int max = 0;
+
+        // For each point in preimage
+        for (int j = 0; j < num; j++) {
+          // Update number of components in preimage
+          if (component[j] > max) max = component[j];
+
+          // Identify component with Cantor polynomial N^2 -> N
+          int identifier = (std::pow(i + component[j], 2) + 3 * i + component[j]) / 2;
+
+          // Update covers
+          cover[preimages[i][j]].push_back(identifier);
+          cover_back[identifier].push_back(preimages[i][j]);
+          cover_fct[identifier] = i;
+          cover_std[identifier] = funcstd[i];
+          cover_color[identifier].second += func_color[preimages[i][j]];
+          cover_color[identifier].first += 1;
+        }
+
+        // Maximal dimension is total number of connected components
+        id += max + 1;
+      }
+    #endif
 
     maximal_dim = id - 1;
     for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++)
@@ -803,24 +852,46 @@ class Cover_complex {
     for (int j = 0; j < n; j++) mindist[j] = std::numeric_limits<double>::max();
 
     // Compute the geodesic distances to subsamples with Dijkstra
-    // #pragma omp parallel for
-    for (int i = 0; i < m; i++) {
-      if (verbose) std::cout << "Computing geodesic distances to seed " << i << "..." << std::endl;
-      int seed = voronoi_subsamples[i];
-      std::vector<double> dmap(n);
-      boost::dijkstra_shortest_paths(
-          one_skeleton, vertices[seed],
-          boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index)));
-
-      for (int j = 0; j < n; j++)
-        if (mindist[j] > dmap[j]) {
-          mindist[j] = dmap[j];
-          if (cover[j].size() == 0)
-            cover[j].push_back(i);
-          else
-            cover[j][0] = i;
-        }
-    }
+    #ifdef GUDHI_USE_TBB
+      if (verbose) std::cout << "Computing geodesic distances (parallelized)..." << std::endl;
+      tbb::mutex coverMutex; tbb::mutex mindistMutex;
+      tbb::parallel_for(0, m, [&](int i){
+        int seed = voronoi_subsamples[i];
+        std::vector<double> dmap(n);
+        boost::dijkstra_shortest_paths(
+            one_skeleton, vertices[seed],
+            boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index)));
+
+        coverMutex.lock(); mindistMutex.lock();
+        for (int j = 0; j < n; j++)
+          if (mindist[j] > dmap[j]) {
+            mindist[j] = dmap[j];
+            if (cover[j].size() == 0)
+              cover[j].push_back(i);
+            else
+              cover[j][0] = i;
+          }
+        coverMutex.unlock(); mindistMutex.unlock();
+      });
+    #else
+      for (int i = 0; i < m; i++) {
+        if (verbose) std::cout << "Computing geodesic distances to seed " << i << "..." << std::endl;
+        int seed = voronoi_subsamples[i];
+        std::vector<double> dmap(n);
+        boost::dijkstra_shortest_paths(
+            one_skeleton, vertices[seed],
+            boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index)));
+
+        for (int j = 0; j < n; j++)
+          if (mindist[j] > dmap[j]) {
+            mindist[j] = dmap[j];
+            if (cover[j].size() == 0)
+              cover[j].push_back(i);
+            else
+              cover[j][0] = i;
+          }
+      }
+    #endif
 
     for (int i = 0; i < n; i++) {
       cover_back[cover[i][0]].push_back(i);
@@ -860,7 +931,7 @@ class Cover_complex {
     while (std::getline(input, line)) {
       std::stringstream stream(line);
       stream >> f;
-      func_color.emplace(i, f);
+      func_color.push_back(f);
       i++;
     }
     color_name = color_file_name;
@@ -873,7 +944,7 @@ class Cover_complex {
            *
            */
   void set_color_from_coordinate(int k = 0) {
-    for (int i = 0; i < n; i++) func_color[i] = point_cloud[i][k];
+    for (int i = 0; i < n; i++) func_color.push_back(point_cloud[i][k]);
     color_name = "coordinate ";
     color_name.append(std::to_string(k));
   }
@@ -885,7 +956,7 @@ class Cover_complex {
    *
    */
   void set_color_from_vector(std::vector<double> color) {
-    for (unsigned int i = 0; i < color.size(); i++) func_color[i] = color[i];
+    for (unsigned int i = 0; i < color.size(); i++) func_color.push_back(color[i]);
   }
 
  public:  // Create a .dot file that can be compiled with neato to produce a .pdf file.
@@ -1039,45 +1110,29 @@ class Cover_complex {
       minf = std::min(minf, it->second);
     }
 
+    // Build filtration
     for (auto const& simplex : simplices) {
-      // Add a simplex and a cone on it
-      std::vector<int> splx = simplex;
-      splx.push_back(-2);
-      st.insert_simplex_and_subfaces(splx);
+      std::vector<int> splx = simplex; splx.push_back(-2);
+      st.insert_simplex_and_subfaces(splx, -3);
     }
 
-    // Build filtration
-    for (auto simplex : st.complex_simplex_range()) {
-      double filta = std::numeric_limits<double>::lowest();
-      double filts = filta;
-      bool ascending = true;
-      for (auto vertex : st.simplex_vertex_range(simplex)) {
-        if (vertex == -2) {
-          ascending = false;
-          continue;
-        }
-        filta = std::max(-2 + (cover_std[vertex] - minf) / (maxf - minf), filta);
-        filts = std::max(2 - (cover_std[vertex] - minf) / (maxf - minf), filts);
-      }
-      if (ascending)
-        st.assign_filtration(simplex, filta);
-      else
-        st.assign_filtration(simplex, filts);
+    for (std::map<int, double>::iterator it = cover_std.begin(); it != cover_std.end(); it++) {
+      int vertex = it->first; float val = it->second;
+      int vert[] = {vertex}; int edge[] = {vertex, -2};
+      st.assign_filtration(st.find(vert), -2 + (val - minf)/(maxf - minf));
+      st.assign_filtration(st.find(edge),  2 - (val - minf)/(maxf - minf));
     }
-    int magic[] = {-2};
-    st.assign_filtration(st.find(magic), -3);
+    st.make_filtration_non_decreasing();
 
     // Compute PD
-    st.initialize_filtration();
-    Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp> pcoh(st);
-    pcoh.init_coefficients(2);
+    Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp> pcoh(st); pcoh.init_coefficients(2);
     pcoh.compute_persistent_cohomology();
 
     // Output PD
     int max_dim = st.dimension();
     for (int i = 0; i < max_dim; i++) {
       std::vector<std::pair<double, double> > bars = pcoh.intervals_in_dimension(i);
-      int num_bars = bars.size();
+      int num_bars = bars.size(); if(i == 0)  num_bars -= 1;
       if(verbose)  std::cout << num_bars << " interval(s) in dimension " << i << ":" << std::endl;
       for (int j = 0; j < num_bars; j++) {
         double birth = bars[j].first;
@@ -1206,8 +1261,7 @@ class Cover_complex {
     }
 
     if (type == "Nerve") {
-      for(auto& simplex : cover)
-        simplices.push_back(simplex.second);
+      for(int i = 0; i < n; i++)  simplices.push_back(cover[i]);
       std::sort(simplices.begin(), simplices.end());
       std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
       simplices.resize(std::distance(simplices.begin(), it));
diff --git a/src/Nerve_GIC/test/test_GIC.cpp b/src/Nerve_GIC/test/test_GIC.cpp
index d633753c..e3067d35 100644
--- a/src/Nerve_GIC/test/test_GIC.cpp
+++ b/src/Nerve_GIC/test/test_GIC.cpp
@@ -39,6 +39,7 @@ BOOST_AUTO_TEST_CASE(check_nerve) {
   N.set_type("Nerve");
   std::string cloud_file_name("data/cloud");
   N.read_point_cloud(cloud_file_name);
+  N.set_color_from_coordinate();
   std::string graph_file_name("data/graph");
   N.set_graph_from_file(graph_file_name);
   std::string cover_file_name("data/cover");
@@ -58,6 +59,7 @@ BOOST_AUTO_TEST_CASE(check_GIC) {
   GIC.set_type("GIC");
   std::string cloud_file_name("data/cloud");
   GIC.read_point_cloud(cloud_file_name);
+  GIC.set_color_from_coordinate();
   std::string graph_file_name("data/graph");
   GIC.set_graph_from_file(graph_file_name);
   std::string cover_file_name("data/cover");
@@ -77,6 +79,7 @@ BOOST_AUTO_TEST_CASE(check_voronoiGIC) {
   GIC.set_type("GIC");
   std::string cloud_file_name("data/cloud");
   GIC.read_point_cloud(cloud_file_name);
+  GIC.set_color_from_coordinate();
   std::string graph_file_name("data/graph");
   GIC.set_graph_from_file(graph_file_name);
   GIC.set_cover_from_Voronoi(Gudhi::Euclidean_distance(), 2);