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

Former-commit-id: 317b90b79f730ea15282d3864ab20d10e148b6e0

2 files changed, 139 insertions, 14 deletions

diff --git a/src/Nerve_GIC/example/simple_GIC.cpp b/src/Nerve_GIC/example/simple_GIC.cpp
index f7a87e85..8b3aecb8 100644
--- a/src/Nerve_GIC/example/simple_GIC.cpp
+++ b/src/Nerve_GIC/example/simple_GIC.cpp
@@ -12,6 +12,7 @@ int main(int argc, char **argv) {
 
   std::string off_file_name(argv[1]);
   double threshold = atof(argv[2]);
+  //std::string func_file_name = argv[3];
   int coord = atoi(argv[3]);
   double resolution = atof(argv[4]);
   double gain = atof(argv[5]);
@@ -30,12 +31,23 @@ int main(int argc, char **argv) {
   GIC.set_graph_from_automatic_rips(100, off_file_name);
   //GIC.set_graph_from_rips(threshold, off_file_name);
   //GIC.set_graph_from_OFF(off_file_name);
+
   GIC.set_function_from_coordinate(coord, off_file_name);
+  //GIC.set_function_from_file(func_file_name);
+
+  GIC.set_color_from_coordinate(coord, off_file_name);
+  //GIC.set_color_from_file(func_file_name);
+
   resolution = GIC.set_automatic_resolution_for_GICMAP();
+
   GIC.set_cover_from_function(resolution,gain,1);
+
   //GIC.find_GIC_simplices();
   //GIC.find_Nerve_simplices();
   GIC.find_GICMAP_simplices_with_functional_minimal_cover();
+
+  GIC.plot_with_KeplerMapper();
+
   Simplex_tree stree; GIC.create_complex(stree);
 
   std::streambuf* streambufffer;
diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h
index 62b5a237..5ee029b7 100644
--- a/src/Nerve_GIC/include/gudhi/GIC.h
+++ b/src/Nerve_GIC/include/gudhi/GIC.h
@@ -48,6 +48,7 @@
 
 #define CONSTANT 10
 #define ETA 0.001
+#define MASK 0
 
 using Simplex_tree = Gudhi::Simplex_tree<>;
 using Filtration_value = Simplex_tree::Filtration_value;
@@ -55,6 +56,7 @@ using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>;
 using Point = std::vector<float>;
 
 std::map<int, double> func;
+std::map<int, double> func_color;
 
 namespace Gudhi {
 
@@ -84,12 +86,15 @@ class Graph_induced_complex {
    std::map<int, std::vector<Cover_t> > cover;
 
  private:
-   int maximal_dim;
+   int maximal_dim; int data_dimension;
 
  private:
    std::map<Cover_t,int> cover_fct;
 
  private:
+  std::map<Cover_t,std::pair<int,double> > cover_color;
+
+ private:
    Simplex_tree<> st;
    std::map<int,std::vector<int> > adjacency_matrix;
 
@@ -176,7 +181,7 @@ class Graph_induced_complex {
    void set_graph_from_rips(const double& threshold, const std::string& off_file_name){
      Points_off_reader<Point> off_reader(off_file_name);
      Rips_complex rips_complex_from_points(off_reader.get_point_cloud(), threshold, Euclidean_distance());
-     rips_complex_from_points.create_complex(st, 1);
+     rips_complex_from_points.create_complex(st, 1); data_dimension = off_reader.get_point_cloud()[0].size();
    }
 
  public: // Automatic tuning of Rips complex.
@@ -185,10 +190,40 @@ class Graph_induced_complex {
      Points_off_reader<Point> off_reader(off_file_name);
      int n = off_reader.get_point_cloud().size();
      int m = floor(n/pow(log(n)/log(CONSTANT),1+ETA)); m = std::min(m,n-1);
-     std::vector<int> samples(m); double delta = 0; int dim = off_reader.get_point_cloud()[0].size();
+     std::vector<int> samples(m); double delta = 0; int dim = off_reader.get_point_cloud()[0].size(); data_dimension = dim;
+
+     std::cout << n << " points in R^" << dim << std::endl;
+     std::cout << "Subsampling " << m << " points" << std::endl;
 
-     std::vector<std::vector<double> > dist; std::vector<double> dumb(n);
-     for(int i = 0; i < n; i++)  dist.push_back(dumb);
+     std::vector<std::vector<double> > dist(n); std::vector<double> dumb(n);
+     for(int i = 0; i < n; i++)  dist[i] = dumb;
+     double d;
+
+     char distances[100]; sprintf(distances,"%s_dist",(char*) off_file_name.c_str());
+     std::ifstream input(distances, std::ios::out | std::ios::binary);
+
+     if(input.good()){
+       std::cout << "Reading distances..." << std::endl;
+       for(int i = 0; i < n; i++){
+         for (int j = 0; j < n; j++){
+           input.read((char*) &d,8); dist[i][j] = d;
+         }
+       }
+       input.close();
+     }
+     else{
+       std::cout << "Computing distances..." << std::endl;
+       input.close(); std::ofstream output(distances, std::ios::out | std::ios::binary);
+       for(int i = 0; i < n; i++){
+         if( (int) floor( 100*(i*1.0+1)/(n*1.0) ) %10 == 0  )
+           std::cout << "\r" << floor( 100*(i*1.0+1)/(n*1.0) ) << "%" << std::flush;
+         for (int j = 0; j < n; j++){
+           double dis = 0; for(int k = 0; k < dim; k++)  dis += pow(off_reader.get_point_cloud()[i][k]-off_reader.get_point_cloud()[j][k],2);
+           dis = std::sqrt(dis); dist[i][j] = dis; output.write((char*) &dis,8);
+         }
+       }
+       output.close(); std::cout << std::endl;
+     }
 
      for(int i = 0; i < n; i++){
        for(int j = i; j < n; j++){
@@ -210,6 +245,7 @@ class Graph_induced_complex {
 
      }
 
+     std::cout << "delta = " << delta << std::endl;
      Rips_complex rips_complex_from_points(off_reader.get_point_cloud(), delta, Euclidean_distance());
      rips_complex_from_points.create_complex(st, 1);
 
@@ -269,6 +305,7 @@ class Graph_induced_complex {
          reso = std::max(reso, std::abs(func[vertices[0]] - func[vertices[1]]));
        }
      }
+     std::cout << "resolution = " << reso << std::endl;
      return reso;
    }
 
@@ -291,7 +328,7 @@ class Graph_induced_complex {
      // Read function values and compute min and max
      std::map<int, double>::iterator it; double maxf, minf; minf = std::numeric_limits<float>::max(); maxf = std::numeric_limits<float>::min();
      for(it = func.begin(); it != func.end(); it++){minf = std::min(minf, it->second); maxf = std::max(maxf, it->second);}
-     int num_pts = func.size(); //std::cout << minf << " " << maxf << std::endl;
+     int num_pts = func.size(); std::cout << "Min function value = " << minf << " and Max function value = " << maxf << std::endl;
 
      // Compute cover of im(f)
      std::vector<std::pair<double,double> > intervals; int res;
@@ -305,7 +342,7 @@ class Graph_induced_complex {
        }
        x = minf + (resolution-1)*incr - alpha; y = maxf;
        std::pair<double, double> interM(x,y); intervals.push_back(interM); res = intervals.size();
-       //for(int i = 0; i < res; i++)  std::cout << intervals[i].first << " " << intervals[i].second << std::endl;
+       for(int i = 0; i < res; i++)  std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" << std::endl;
      }
      else{
        double x = minf; double y = x + resolution;
@@ -315,6 +352,7 @@ class Graph_induced_complex {
          y = x + resolution;
        }
        std::pair<double, double> interM(x,maxf); intervals.push_back(interM); res = intervals.size();
+       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
@@ -382,27 +420,102 @@ class Graph_induced_complex {
        }
 
        // Compute the connected components with DFS
-       std::map<int,bool> visit; //std::cout << i << std::endl;
+       std::map<int,bool> visit; std::cout << "Preimage of interval " << i << std::endl;
        for(std::map<int, std::vector<int> >::iterator it = prop.begin(); it != prop.end(); it++)
          visit.insert(std::pair<int,bool>(it->first, false));
        if (!(prop.empty())){
          for(std::map<int, std::vector<int> >::iterator it = prop.begin(); it != prop.end(); it++){
            if (  !(visit[it->first])  ){
              std::vector<int> cc; cc.clear();
-             dfs(prop,it->first,cc,visit); int cci = cc.size(); //std::cout << cci << " ";
-             for(int j = 0; j < cci; j++)  cover[cc[j]].push_back(id);
-             cover_fct[id] = i;
+             dfs(prop,it->first,cc,visit); int cci = cc.size(); std::cout << "one CC with " << cci << " points, ";
+             double average_col = 0;
+             for(int j = 0; j < cci; j++){cover[cc[j]].push_back(id); average_col += func_color[cc[j]]/cci;}
+             cover_fct[id] = i; cover_color[id] = std::pair<int,double>(cci,average_col);
              id++;
            }
          }
        }
-       //std::cout << std::endl;
+       std::cout << std::endl;
      }
 
      maximal_dim = id;
 
    }
 
+ // *******************************************************************************************************************
+ // Visualization.
+ // *******************************************************************************************************************
+
+ public: // Set color from file.
+  void set_color_from_file(const std::string& color_file_name){
+    int vertex_id = 0; std::ifstream input(color_file_name); std::string line; double f;
+    while(std::getline(input,line)){
+      std::stringstream stream(line); stream >> f;
+      func_color.insert(std::pair<int,double>(vertex_id, f)); vertex_id++;
+    }
+  }
+
+ public: // Set color from kth coordinate
+  void set_color_from_coordinate(const int& k, const std::string& off_file_name){
+    Points_off_reader<Point> off_reader(off_file_name);
+    int numpts = off_reader.get_point_cloud().size();
+    for(int i = 0; i < numpts; i++)  func_color.insert(std::pair<int,double>(i,off_reader.get_point_cloud()[i][k]));
+  }
+
+ public: // Set color from vector.
+  void set_color_from_vector(const std::vector<double>& color){
+    for(int i = 0; i < color.size(); i++)  func_color.insert(std::pair<int,double>(i, color[i]));
+  }
+
+ public: // Create a .dot file that can be compiled with neato to produce a .pdf file
+  void plot_with_neato(){
+    char mapp[11] = "mapper.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++){
+      maxv = std::max(maxv, iit->second.second);  minv = std::min(minv, iit->second.second);
+    }
+    int k = 0; std::vector<int> nodes; nodes.clear();
+    for (std::map<Cover_t,std::pair<int,double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++){
+      if(iit->second.first > MASK){
+        nodes.push_back(iit->first);
+        graphic << iit->first << "[shape=circle fontcolor=black color=black label=\"" \
+          << iit->first << ":" << iit->second.first << "\" style=filled fillcolor=\"" \
+          << (1-(maxv-iit->second.second)/(maxv-minv))*0.6 << ", 1, 1\"]" << std::endl;
+        k++;
+      }
+    }
+    int ke = 0; int num_simplices = simplices.size();
+    for (int i = 0; i < num_simplices; i++)
+      if (simplices[i].size() == 2)
+        if (cover_color[simplices[i][0]].first > MASK && cover_color[simplices[i][1]].first > MASK){
+          graphic << "  " << simplices[i][0] << " -- " << simplices[i][1] << " [weight=15];" << std::endl; ke++;}
+    graphic << "}"; graphic.close();
+  }
+
+ public: // Create a .txt file that can be compiled with KeplerMapper to produce a .html file
+  void plot_with_KeplerMapper(){
+
+    int num_simplices = simplices.size(); int num_edges = 0;
+    char mapp[11] = "mapper.txt";  std::ofstream graphic(mapp);
+    for (int i = 0; i < num_simplices; i++)
+      if (simplices[i].size() == 2)
+        if (cover_color[simplices[i][0]].first > MASK && cover_color[simplices[i][1]].first > MASK)
+          num_edges++;
+
+    graphic << "Cloud" << std::endl;
+    graphic << "Function" << std::endl;
+    graphic << 0 << " " << 0 << std::endl;
+    graphic << cover_color.size() << " " << num_edges << std::endl;
+
+    for (std::map<Cover_t,std::pair<int,double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++)
+      graphic << iit->first << " " << iit->second.second << " " << iit->second.first << std::endl;
+
+    for (int i = 0; i < num_simplices; i++)
+      if (simplices[i].size() == 2)
+        if (cover_color[simplices[i][0]].first > MASK && cover_color[simplices[i][1]].first > MASK)
+          graphic << simplices[i][0] << " " << simplices[i][1] << std::endl;
+    graphic.close();
+  }
 
  // *******************************************************************************************************************
  // *******************************************************************************************************************
@@ -490,10 +603,10 @@ class Graph_induced_complex {
      for(std::map<int,std::vector<Cover_t> >::iterator it = cover.begin(); it != cover.end(); it++){
 
        int vid = it->first; std::vector<int> neighbors = adjacency_matrix[vid]; int num_neighb = neighbors.size();
-       //std::cout << vid << " " << num_neighb << std::endl;
 
        // Find cover of current point (vid).
        if(cover[vid].size() == 2)  v1 = std::min(cover[vid][0],cover[vid][1]); else  v1 = cover[vid][0];
+       std::vector<int> node(1); node[0] = v1; simplices.push_back(node);
 
        // Loop on neighbors.
        for(int i = 0; i < num_neighb; i++){
@@ -505,7 +618,7 @@ class Graph_induced_complex {
 
          // If neighbor is in next interval, add edge.
          if(cover_fct[v2] == cover_fct[v1] + 1){
-           std::vector<int> edge(2); edge[0] = v1; edge[1] = v2; //std::cout << v1 << " " << v2 << std::endl;
+           std::vector<int> edge(2); edge[0] = v1; edge[1] = v2;
            simplices.push_back(edge);
          }
        }