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

Former-commit-id: 1624f24a6f87ac092e1761c9778741102fb37a3c

7 files changed, 132 insertions, 105 deletions

diff --git a/src/Nerve_GIC/doc/Intro_graph_induced_complex.h b/src/Nerve_GIC/doc/Intro_graph_induced_complex.h
index 7120ada3..8d807f4b 100644
--- a/src/Nerve_GIC/doc/Intro_graph_induced_complex.h
+++ b/src/Nerve_GIC/doc/Intro_graph_induced_complex.h
@@ -164,7 +164,7 @@ namespace graph_induced_complex {
  * where the graph G comes from a Rips complex with optimal threshold,
  * and the cover C comes from the preimages of intervals covering the first coordinate,
  * with optimal resolution and gain. Note that optimal threshold, resolution and gain
- * also exist for the Nerve of this cover.
+ * can be computed as well for the Nerve.
  *
  * \include Nerve_GIC/MapperDeltaCoord.cpp
  *
diff --git a/src/Nerve_GIC/example/GIC.cpp b/src/Nerve_GIC/example/GIC.cpp
index aba0927d..6e5c5ca7 100644
--- a/src/Nerve_GIC/example/GIC.cpp
+++ b/src/Nerve_GIC/example/GIC.cpp
@@ -30,7 +30,7 @@ void usage(int nbArgs, char * const progName) {
 }
 
 int main(int argc, char **argv) {
-  if ((argc != 6) && (argc != 7)) usage(argc, (argv[0] - 1));
+  if ((argc != 6) && (argc != 7)) usage(argc, argv[0]);
 
   using Point = std::vector<float>;
 
@@ -58,33 +58,30 @@ int main(int argc, char **argv) {
 
     GIC.set_graph_from_rips(threshold, Gudhi::Euclidean_distance());
 
-    GIC.set_resolution_double(resolution); GIC.set_gain(gain);
-    GIC.set_cover_from_function(1);
+    GIC.set_resolution_with_interval_length(resolution); GIC.set_gain(gain);
+    GIC.set_cover_from_function();
 
     GIC.find_GIC_simplices();
 
-    GIC.plot_txt();
+    GIC.plot_TXT_for_KeplerMapper();
 
     Gudhi::graph_induced_complex::Simplex_tree stree; GIC.create_complex(stree);
 
-    std::streambuf* streambufffer = std::cout.rdbuf();
-    std::ostream output_stream(streambufffer);
-
     // ----------------------------------------------------------------------------
     // Display information about the graph induced complex
     // ----------------------------------------------------------------------------
 
     if(verb){
-      output_stream << "Graph induced complex is of dimension " << stree.dimension() <<
+      std::cout << "Graph induced complex is of dimension " << stree.dimension() <<
                    " - " << stree.num_simplices() << " simplices - " <<
                    stree.num_vertices() << " vertices." << std::endl;
 
-      output_stream << "Iterator on graph induced complex simplices" << std::endl;
+      std::cout << "Iterator on graph induced complex simplices" << std::endl;
       for (auto f_simplex : stree.filtration_simplex_range()) {
         for (auto vertex : stree.simplex_vertex_range(f_simplex)) {
-          output_stream << vertex << " ";
+          std::cout << vertex << " ";
         }
-        output_stream << std::endl;
+        std::cout << std::endl;
       }
     }
   }
diff --git a/src/Nerve_GIC/example/GICvoronoi.cpp b/src/Nerve_GIC/example/GICvoronoi.cpp
index 013d6290..2c4f5acf 100644
--- a/src/Nerve_GIC/example/GICvoronoi.cpp
+++ b/src/Nerve_GIC/example/GICvoronoi.cpp
@@ -30,7 +30,7 @@ void usage(int nbArgs, char * const progName) {
 }
 
 int main(int argc, char **argv) {
-  if ((argc != 3) && (argc != 4)) usage(argc, (argv[0] - 1));
+  if ((argc != 3) && (argc != 4)) usage(argc, argv[0]);
 
   using Point = std::vector<float>;
 
@@ -52,7 +52,7 @@ int main(int argc, char **argv) {
 
     GIC.set_color_from_coordinate();
 
-    GIC.set_graph_from_OFF(off_file_name);
+    GIC.set_graph_from_OFF();
 
     GIC.set_cover_from_Voronoi(Gudhi::Euclidean_distance(),m);
 
@@ -62,24 +62,21 @@ int main(int argc, char **argv) {
 
     Gudhi::graph_induced_complex::Simplex_tree stree; GIC.create_complex(stree);
 
-    std::streambuf* streambufffer = std::cout.rdbuf();
-    std::ostream output_stream(streambufffer);
-
     // ----------------------------------------------------------------------------
     // Display information about the graph induced complex
     // ----------------------------------------------------------------------------
 
     if(verb){
-      output_stream << "Graph induced complex is of dimension " << stree.dimension() <<
+      std::cout << "Graph induced complex is of dimension " << stree.dimension() <<
                      " - " << stree.num_simplices() << " simplices - " <<
                      stree.num_vertices() << " vertices." << std::endl;
 
-      output_stream << "Iterator on graph induced complex simplices" << std::endl;
+      std::cout << "Iterator on graph induced complex simplices" << std::endl;
       for (auto f_simplex : stree.filtration_simplex_range()) {
         for (auto vertex : stree.simplex_vertex_range(f_simplex)) {
-          output_stream << vertex << " ";
+          std::cout << vertex << " ";
         }
-        output_stream << std::endl;
+        std::cout << std::endl;
       }
     }
   }
diff --git a/src/Nerve_GIC/example/MapperDeltaCoord.cpp b/src/Nerve_GIC/example/MapperDeltaCoord.cpp
index 9927f762..382649e8 100644
--- a/src/Nerve_GIC/example/MapperDeltaCoord.cpp
+++ b/src/Nerve_GIC/example/MapperDeltaCoord.cpp
@@ -30,7 +30,7 @@ void usage(int nbArgs, char * const progName) {
 }
 
 int main(int argc, char **argv) {
-  if ((argc != 3) && (argc != 4)) usage(argc, (argv[0] - 1));
+  if ((argc != 3) && (argc != 4)) usage(argc, argv[0]);
 
   using Point = std::vector<float>;
 
@@ -56,32 +56,29 @@ int main(int argc, char **argv) {
     GIC.set_graph_from_automatic_rips(Gudhi::Euclidean_distance());
 
     GIC.set_automatic_resolution_for_GICMAP(); GIC.set_gain();
-    GIC.set_cover_from_function(1);
+    GIC.set_cover_from_function();
 
     GIC.find_GICMAP_simplices_with_functional_minimal_cover();
 
-    GIC.plot_pdf();
+    GIC.plot_DOT_for_neato();
 
     Gudhi::graph_induced_complex::Simplex_tree stree; GIC.create_complex(stree);
 
-    std::streambuf* streambufffer = std::cout.rdbuf();
-    std::ostream output_stream(streambufffer);
-
     // ------------------------------------------
     // Display information about the Mapper Delta
     // ------------------------------------------
 
     if(verb){
-      output_stream << "Mapper Delta is of dimension " << stree.dimension() <<
+      std::cout << "Mapper Delta is of dimension " << stree.dimension() <<
                      " - " << stree.num_simplices() << " simplices - " <<
                      stree.num_vertices() << " vertices." << std::endl;
 
-      output_stream << "Iterator on Mapper Delta simplices" << std::endl;
+      std::cout << "Iterator on Mapper Delta simplices" << std::endl;
       for (auto f_simplex : stree.filtration_simplex_range()) {
         for (auto vertex : stree.simplex_vertex_range(f_simplex)) {
-          output_stream << vertex << " ";
+          std::cout << vertex << " ";
         }
-        output_stream << std::endl;
+        std::cout << std::endl;
       }
     }
   }
diff --git a/src/Nerve_GIC/example/MapperDeltaFunc.cpp b/src/Nerve_GIC/example/MapperDeltaFunc.cpp
index f3d57d2f..586b733c 100644
--- a/src/Nerve_GIC/example/MapperDeltaFunc.cpp
+++ b/src/Nerve_GIC/example/MapperDeltaFunc.cpp
@@ -30,7 +30,7 @@ void usage(int nbArgs, char * const progName) {
 }
 
 int main(int argc, char **argv) {
-  if ((argc != 3) && (argc != 4)) usage(argc, (argv[0] - 1));
+  if ((argc != 3) && (argc != 4)) usage(argc, argv[0]);
 
   using Point = std::vector<float>;
 
@@ -56,32 +56,29 @@ int main(int argc, char **argv) {
     GIC.set_graph_from_automatic_rips(Gudhi::Euclidean_distance());
 
     GIC.set_automatic_resolution_for_GICMAP(); GIC.set_gain();
-    GIC.set_cover_from_function(1);
+    GIC.set_cover_from_function();
 
     GIC.find_GICMAP_simplices_with_functional_minimal_cover();
 
-    GIC.plot_pdf();
+    GIC.plot_DOT_for_neato();
 
     Gudhi::graph_induced_complex::Simplex_tree stree; GIC.create_complex(stree);
 
-    std::streambuf* streambufffer = std::cout.rdbuf();
-    std::ostream output_stream(streambufffer);
-
     // ------------------------------------------
     // Display information about the Mapper Delta
     // ------------------------------------------
 
     if(verb){
-      output_stream << "Mapper Delta is of dimension " << stree.dimension() <<
+      std::cout << "Mapper Delta is of dimension " << stree.dimension() <<
                      " - " << stree.num_simplices() << " simplices - " <<
                      stree.num_vertices() << " vertices." << std::endl;
 
-      output_stream << "Iterator on Mapper Delta simplices" << std::endl;
+      std::cout << "Iterator on Mapper Delta simplices" << std::endl;
       for (auto f_simplex : stree.filtration_simplex_range()) {
         for (auto vertex : stree.simplex_vertex_range(f_simplex)) {
-          output_stream << vertex << " ";
+          std::cout << vertex << " ";
         }
-        output_stream << std::endl;
+        std::cout << std::endl;
       }
     }
   }
diff --git a/src/Nerve_GIC/example/Nerve.cpp b/src/Nerve_GIC/example/Nerve.cpp
index 10c5ba3b..e47982c4 100644
--- a/src/Nerve_GIC/example/Nerve.cpp
+++ b/src/Nerve_GIC/example/Nerve.cpp
@@ -30,7 +30,7 @@ void usage(int nbArgs, char * const progName) {
 }
 
 int main(int argc, char **argv) {
-  if ((argc != 5) && (argc != 6)) usage(argc, (argv[0] - 1));
+  if ((argc != 5) && (argc != 6)) usage(argc, argv[0]);
 
   using Point = std::vector<float>;
 
@@ -55,35 +55,32 @@ int main(int argc, char **argv) {
     GIC.set_color_from_coordinate(coord);
     GIC.set_function_from_coordinate(coord);
 
-    GIC.set_graph_from_OFF(off_file_name);
+    GIC.set_graph_from_OFF();
 
-    GIC.set_resolution_int(resolution); GIC.set_gain(gain);
-    GIC.set_cover_from_function(0);
+    GIC.set_resolution_with_interval_number(resolution); GIC.set_gain(gain);
+    GIC.set_cover_from_function();
 
     GIC.find_Nerve_simplices();
 
-    GIC.plot_txt();
+    GIC.plot_TXT_for_KeplerMapper();
 
     Gudhi::graph_induced_complex::Simplex_tree stree; GIC.create_complex(stree);
 
-    std::streambuf* streambufffer = std::cout.rdbuf();
-    std::ostream output_stream(streambufffer);
-
     // ----------------------------------------------------------------------------
     // Display information about the graph induced complex
     // ----------------------------------------------------------------------------
 
     if(verb){
-      output_stream << "Nerve is of dimension " << stree.dimension() <<
+      std::cout << "Nerve is of dimension " << stree.dimension() <<
                      " - " << stree.num_simplices() << " simplices - " <<
                      stree.num_vertices() << " vertices." << std::endl;
 
-      output_stream << "Iterator on Nerve simplices" << std::endl;
+      std::cout << "Iterator on Nerve simplices" << std::endl;
       for (auto f_simplex : stree.filtration_simplex_range()) {
         for (auto vertex : stree.simplex_vertex_range(f_simplex)) {
-          output_stream << vertex << " ";
+          std::cout << vertex << " ";
         }
-        output_stream << std::endl;
+        std::cout << std::endl;
       }
     }
   }
diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h
index 42225c47..f00856ed 100644
--- a/src/Nerve_GIC/include/gudhi/GIC.h
+++ b/src/Nerve_GIC/include/gudhi/GIC.h
@@ -77,6 +77,7 @@ class Graph_induced_complex {
    //Graph_induced_complex(std::map<int, double> fun){func = fun;}
    bool verbose; // whether to display information.
    std::vector<Point> point_cloud;
+   std::vector<std::vector<int> > one_skeleton;
    typedef int Cover_t; // elements of cover C are indexed by integers.
    std::vector<std::vector<Cover_t> > simplices;
    std::map<int, std::vector<Cover_t> > cover;
@@ -88,9 +89,9 @@ class Graph_induced_complex {
    Simplex_tree st;
    std::map<int,std::vector<int> > adjacency_matrix;
    std::vector<std::vector<double> > distances;
-   int resolution_int;
-   double resolution_double;
-   double gain;
+   int resolution_int = -1;
+   double resolution_double = -1;
+   double gain = -1;
    double rate_constant; // Constant in the subsampling.
    double rate_power; // Power in the subsampling.
    int mask; // Ignore nodes containing less than mask points.
@@ -152,23 +153,62 @@ class Graph_induced_complex {
    }
 
  public:
+   /** \brief Specifies whether the program should display information or not.
+    *
+    * @param[in] verb boolean (true = display info, false = do not display info).
+    *
+    */
    void set_verbose(bool verb = 0){verbose = verb;}
  public:
+   /** \brief Sets the constants used to subsample the data set. These constants are
+    * explained in "Statistical Analysis and Parameter Selection for the Mapper".
+    *
+    * @param[in] constant double.
+    * @param[in] power double.
+    *
+    */
    void set_subsampling(double constant = 10, double power = 0.001){rate_constant = constant; rate_power = power;}
  public:
+   /** \brief Sets the mask, which is a threshold integer such that nodes in the complex that contain less data points
+    * than this threshold are not displayed.
+    *
+    * @param[in] nodemask integer.
+    *
+    */
    void set_mask(int nodemask = 0){mask = nodemask;}
 
  public:
+   /** \brief Reads and stores the input point cloud.
+    *
+    * @param[in] off_file_name name of the input .OFF file.
+    *
+    */
    bool read_point_cloud(std::string off_file_name){
-     Gudhi::Points_off_reader<Point> off_reader = Points_off_reader<Point>(off_file_name);
-     bool check = off_reader.is_valid();
-     if(check){
-       point_cloud = off_reader.get_point_cloud();
-       point_cloud_name = off_file_name;
-       n = point_cloud.size();
-       data_dimension = point_cloud[0].size();
-       return check;
-     } else  return check;
+     point_cloud_name = off_file_name;
+     int numedges, numfaces, i, num; std::vector<int> edge(2);
+     std::vector<int> simplex;
+     std::ifstream input(off_file_name); std::string line; getline(input, line);
+     input >> n; input >> numfaces; input >> numedges; getline(input, line);
+
+     i = 0; while(i < n){
+       getline(input, line);
+       std::vector<double> point; std::istringstream iss(line);
+       point.assign(std::istream_iterator<double>(iss), std::istream_iterator<double>());
+       point_cloud.push_back(Point(point.begin(),point.end())); i++;
+     } data_dimension = point_cloud[0].size();
+
+     i = 0;  while(i < numfaces){
+       simplex.clear(); input >> num;
+       for(int j = 0; j < num; j++){int k; input >> k; simplex.push_back(k);}
+       for(int j = 0; j < num; j++){
+         for(int k = j+1; k < num; k++){
+           edge[0] = simplex[j]; edge[1] = simplex[k]; one_skeleton.push_back(edge);
+         }
+       }
+       i++;
+     }
+
+     return input.is_open();
    }
 
  // *******************************************************************************************************************
@@ -197,29 +237,17 @@ class Graph_induced_complex {
    }
 
  public: // Set graph from OFF file.
-   /** \brief Creates the graph G from the triangulation given by an .OFF file.
-    *
-    * @param[in] off_file_name name of the input .OFF file.
+   /** \brief Creates the graph G from the triangulation given by the input .OFF file.
     *
     */
-   void set_graph_from_OFF(std::string off_file_name){
-     int numedges, numfaces, i; std::vector<int> edge(2); double x; int num; std::vector<int> simplex;
-     std::ifstream input(off_file_name); std::string line; getline(input, line);
-     input >> n; input >> numfaces; input >> numedges;
-     i = 0;  while(i < n){input >> x; input >> x; input >> x; i++;}
-     i = 0;  while(i < numfaces){
-       simplex.clear(); input >> num;
-       for(int j = 0; j < num; j++){int k; input >> k; simplex.push_back(k);}
-       for(int j = 0; j < num; j++){
-         for(int k = j+1; k < num; k++){
-           edge[0] = simplex[j]; edge[1] = simplex[k];
-           st.insert_simplex_and_subfaces(edge);
-         }
-       }
-       i++;
-     }
+   void set_graph_from_OFF(){
 
-     fill_adjacency_matrix_from_st();
+     int num_edges = one_skeleton.size();
+     if(num_edges > 0){
+       for(int i = 0; i < num_edges; i++)  st.insert_simplex_and_subfaces(one_skeleton[i]);
+       fill_adjacency_matrix_from_st();
+     }
+     else  std::cout << "No triangulation read in OFF file!" << std::endl;
 
    }
 
@@ -378,13 +406,13 @@ class Graph_induced_complex {
     * @param[in] reso length of intervals.
     *
     */
-   void set_resolution_double(double reso){resolution_double = reso;}
+   void set_resolution_with_interval_length(double reso){resolution_double = reso;}
    /** \brief Sets a number of intervals from a value stored in memory.
     *
     * @param[in] reso number of intervals.
     *
     */
-   void set_resolution_int(int reso){resolution_int = reso;}
+   void set_resolution_with_interval_number(int reso){resolution_int = reso;}
    /** \brief Sets a gain from a value stored in memory (default value 0.3).
     *
     * @param[in] g gain.
@@ -411,10 +439,11 @@ class Graph_induced_complex {
  public: // Set cover with preimages of function.
    /** \brief Creates a cover C from the preimages of the function f.
     *
-    * @param[in] token boolean specifying whether we use the length or the number of intervals for the cover of im(f).
-    *
     */
-   void set_cover_from_function(bool token){
+   void set_cover_from_function(){
+
+     if(resolution_double == -1 && resolution_int == -1){ std::cout << "Number and/or length of intervals not specified" << std::endl; return; }
+     if(gain == -1){ std::cout << "Gain not specified" << std::endl; return; }
 
      // Read function values and compute min and max
      std::map<int, double>::iterator it;
@@ -425,7 +454,7 @@ class Graph_induced_complex {
      // Compute cover of im(f)
      std::vector<std::pair<double,double> > intervals; int res;
 
-     if(!token){ // Case we use an integer for the number of intervals.
+     if(resolution_double == -1){ // Case we use an integer for the number of intervals.
        double incr = (maxf-minf)/resolution_int; double x = minf; double alpha = (incr*gain)/(2-2*gain);
        double y = minf + incr + alpha; std::pair<double, double> interm(x,y); intervals.push_back(interm);
        for(int i = 1; i < resolution_int-1; i++){
@@ -439,16 +468,30 @@ class Graph_induced_complex {
          for(int i = 0; i < res; i++)  std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" << std::endl;
      }
 
-     else{ // Case we use a double for the length of the intervals.
-       double x = minf; double y = x + resolution_double;
-       while(y <= maxf && maxf - (y-gain*resolution_double) >= resolution_double){
-         std::pair<double, double> inter(x,y); intervals.push_back(inter);
-         x = y - gain*resolution_double;
-         y = x + resolution_double;
+     else{
+       if(resolution_int == -1){ // Case we use a double for the length of the intervals.
+         double x = minf; double y = x + resolution_double;
+         while(y <= maxf && maxf - (y-gain*resolution_double) >= resolution_double){
+           std::pair<double, double> inter(x,y); intervals.push_back(inter);
+           x = y - gain*resolution_double;
+           y = x + resolution_double;
+         }
+         std::pair<double, double> interM(x,maxf); intervals.push_back(interM); res = intervals.size();
+         if(verbose)
+           for(int i = 0; i < res; i++)  std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" << std::endl;
+       }
+
+       else{ // Case we use an integer and a double for the length of the intervals.
+         double x = minf; double y = x + resolution_double; int count = 0;
+         while(count < resolution_int && y <= maxf && maxf - (y-gain*resolution_double) >= resolution_double){
+             std::pair<double, double> inter(x,y); intervals.push_back(inter);  count++;
+             x = y - gain*resolution_double;
+             y = x + resolution_double;
+           }
+           res = intervals.size();
+           if(verbose)
+             for(int i = 0; i < res; i++)  std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" << std::endl;
        }
-       std::pair<double, double> interM(x,maxf); intervals.push_back(interM); res = intervals.size();
-       if(verbose)
-         for(int i = 0; i < res; i++)  std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" << std::endl;
      }
 
      // Sort points according to function values
@@ -624,7 +667,6 @@ class Graph_induced_complex {
    /** \brief Computes the function used to color the nodes of the simplicial complex from the k-th coordinate.
     *
     * @param[in] k coordinate to use (start at 0).
-    * @param[in] off_file_name name of the input .OFF file.
     *
     */
    void set_color_from_coordinate(int k = 0){
@@ -643,10 +685,10 @@ class Graph_induced_complex {
    }
 
  public: // Create a .dot file that can be compiled with neato to produce a .pdf file.
-   /** \brief Creates a .dot file for neato once the simplicial complex is computed to get a nice visualization
+   /** \brief Creates a .dot file for neato (part of the graphviz package) once the simplicial complex is computed to get a visualization
     * of its 1-skeleton in a .pdf file.
     */
-   void plot_pdf(){
+   void plot_DOT_for_neato(){
      char mapp[11] = "SC.dot";   std::ofstream graphic(mapp); graphic << "graph Mapper {" << std::endl;
      double maxv, minv; maxv = std::numeric_limits<double>::min(); minv = std::numeric_limits<double>::max();
      for (std::map<Cover_t,std::pair<int,double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++){
@@ -672,10 +714,10 @@ class Graph_induced_complex {
    }
 
  public: // Create a .txt file that can be compiled with KeplerMapper to produce a .html file.
-   /** \brief Creates a .txt file for KeplerMapper once the simplicial complex is computed to get a nice visualization
+   /** \brief Creates a .txt file for KeplerMapper once the simplicial complex is computed to get a visualization
     * of its 1-skeleton in browser.
     */
-   void plot_txt(){
+   void plot_TXT_for_KeplerMapper(){
 
      int num_simplices = simplices.size(); int num_edges = 0;
      char mapp[11] = "SC.txt";  std::ofstream graphic(mapp);