From b3654bd1ac7d38aaac0550de063158dbdf96522f Mon Sep 17 00:00:00 2001
From: mcarrier <mcarrier@636b058d-ea47-450e-bf9e-a15bfbe3eedb>
Date: Tue, 24 Jan 2017 17:06:57 +0000
Subject: added headers and examples for Nerve GIC

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

Former-commit-id: 587a3e6748bb2b7313dbd1bd4e81c4a3b01a5efc
---
 src/Nerve_GIC/example/bottleneckBootstrap.cpp | 274 ++++++++++
 src/Nerve_GIC/example/density_smoother.cpp    | 155 ++++++
 src/Nerve_GIC/example/graph_off.cpp           |  49 ++
 src/Nerve_GIC/example/km.py                   | 390 ++++++++++++++
 src/Nerve_GIC/example/mapper.cpp              | 179 +++++++
 src/Nerve_GIC/example/mapper_parameter.cpp    | 165 ++++++
 src/Nerve_GIC/example/quantile.cpp            |  48 ++
 src/Nerve_GIC/example/visu.py                 |  43 ++
 src/Nerve_GIC/include/gudhi/mapper.h          | 705 ++++++++++++++++++++++++++
 src/Nerve_GIC/include/gudhi/utils.h           |  86 ++++
 10 files changed, 2094 insertions(+)
 create mode 100644 src/Nerve_GIC/example/bottleneckBootstrap.cpp
 create mode 100644 src/Nerve_GIC/example/density_smoother.cpp
 create mode 100644 src/Nerve_GIC/example/graph_off.cpp
 create mode 100755 src/Nerve_GIC/example/km.py
 create mode 100755 src/Nerve_GIC/example/mapper.cpp
 create mode 100644 src/Nerve_GIC/example/mapper_parameter.cpp
 create mode 100644 src/Nerve_GIC/example/quantile.cpp
 create mode 100755 src/Nerve_GIC/example/visu.py
 create mode 100644 src/Nerve_GIC/include/gudhi/mapper.h
 create mode 100644 src/Nerve_GIC/include/gudhi/utils.h

diff --git a/src/Nerve_GIC/example/bottleneckBootstrap.cpp b/src/Nerve_GIC/example/bottleneckBootstrap.cpp
new file mode 100644
index 00000000..eac94e5e
--- /dev/null
+++ b/src/Nerve_GIC/example/bottleneckBootstrap.cpp
@@ -0,0 +1,274 @@
+ #include "mapper.h"
+#define ETA 0.001
+#define SHIFT 1
+#define DELTA 0
+#define CONSTANT 10
+
+double GetUniform(){
+  static default_random_engine re;
+  static uniform_real_distribution<double> Dist(0,1);
+  return Dist(re);
+}
+
+void SampleWithReplacement(int populationSize, int sampleSize, vector<int> & samples){
+  int m = 0; double u;
+  while (m < sampleSize){
+    u = GetUniform();
+    samples[m] = floor(u*populationSize);
+    m++;
+  }
+}
+
+void SampleWithoutReplacement(int populationSize, int sampleSize, vector<int> & samples){
+  int& n = sampleSize; int& N = populationSize;
+  int t = 0; int m = 0; double u;
+  while (m < n){
+    u = GetUniform();
+    if ( (N - t)*u >= n - m )
+      t++;
+    else{samples[m] = t; t++; m++;}
+  }
+}
+
+void compute_estimator(Cloud & C, char* const & cloud, const bool & VNE, const int & Nsubs, \
+                       const double & gain){
+
+    int n = C.size();
+    double** dist = new double*[n];
+    for(int i = 0; i < n; i++)  dist[i] = new double[n];
+    double d;
+    Cover I; AdjacencyMatrix G; MapperElements M;
+    double maxf = C[0].func; double minf = C[0].func;
+    for(int i = 1; i < n; i++){ maxf = max(maxf,C[i].func); minf = min(minf,C[i].func);}
+
+    char name1[100];
+    if(VNE)  sprintf(name1, "%s_VNdist", cloud);
+    else  sprintf(name1, "%s_dist", cloud);
+    char* const name2 = name1;
+    ifstream input(name2, ios::out | ios::binary);
+
+    if(input.good()){
+      cout << "  Reading distances for parameter selection..." << 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{
+      cout << "  Computing distances for parameter selection..." << endl; vector<double> V;
+      input.close(); ofstream output(name2, ios::out | ios::binary);
+      if(VNE){
+        int dim = C[0].coord.size();
+        for(int i = 0; i < dim; i++){
+          double meani = 0;
+          for (int j = 0; j < n; j++)  meani += C[j].coord[i]/n;
+          double vari = 0;
+          for (int j = 0; j < n; j++)  vari += pow(C[j].coord[i]-meani,2)/n;
+          V.push_back(vari);
+        }
+      }
+
+      for(int i = 0; i < n; i++){
+        if( (int) floor( 100*((double) i)/((double) n)+1 ) %10 == 0  ){
+          cout << "\r" << floor( 100*((double) i)/((double) n) +1) << "%" << flush;}
+        for (int j = 0; j < n; j++){
+          if(!VNE){d = C[i].EuclideanDistance(C[j]); dist[i][j] = d;}
+          else{d = C[i].VarianceNormalizedEuclideanDistance(C[j],V); dist[i][j] = d;}
+          output.write((char*) &d,8);
+        }
+      }
+      output.close();
+    }
+
+    int m = floor(n/pow(log(n)/log(CONSTANT),1+ETA));
+    double lip = 0; double delta = 0; vector<int> samples(m);
+    #pragma omp parallel for
+    for (int i = 0; i < Nsubs; i++){
+      SampleWithoutReplacement(n,m,samples);
+      double hausdorff_dist = 0;
+      for (int j = 0; j < n; j++){
+        Point p = C[j]; double mj = dist[j][samples[0]]; double Mi = abs(p.func-C[0].func)/dist[j][0];
+        for (int k = 1; k < m; k++)  mj = min(mj, dist[j][samples[k]]);
+        if(j < n-1)
+          for (int k = j+1; k < n; k++)  Mi = max(Mi, abs(p.func-C[k].func)/dist[j][k]);
+        hausdorff_dist = max(hausdorff_dist, mj); lip = max(lip,Mi);
+      }
+      delta += hausdorff_dist;
+    }
+    delta /= Nsubs; double resolution; if(DELTA)  resolution = lip*delta; else  resolution = lip*delta/gain;
+    //cout << delta << " " << lip << " " << resolution << endl;
+    cout << "  Done." << endl;
+
+    G = build_neighborhood_graph_from_scratch_with_parameter(C,delta,name2,VNE);
+
+    sort(C.begin(),C.end());
+    I = Cover(C.begin()->func, (C.end()-1)->func, resolution, gain, SHIFT, 1);
+    I.proper_value();
+
+    map<int,double> colors; colors.clear(); map<int,int> num; num.clear();
+    map<Point,vector<int> > clusters; vector<int> dum; dum.clear();
+    for(int i = 0; i < C.size(); i++)  clusters.insert(pair<Point,vector<int> >(C[i],dum)); vector<double> col(G.size());
+
+    M = MapperElts(G,I,clusters,col);
+    SparseGraph MG;
+    if (DELTA)  MG = build_mapperDelta_graph(M,G,clusters,colors,num); else  MG = build_mapper_graph(M,colors,num);
+
+    char mappg[100]; sprintf(mappg,"%s_mapper.txt",cloud);
+    ofstream graphicg(mappg);
+
+    double maxv, minv; maxv = numeric_limits<double>::min(); minv = numeric_limits<double>::max();
+    for (map<int,double>::iterator iit = colors.begin(); iit != colors.end(); iit++){
+      if(iit->second > maxv)  maxv = iit->second;
+      if(minv > iit->second)  minv = iit->second;
+    }
+
+    int k = 0; int ke = 0;
+
+    for (SparseGraph::iterator it = MG.begin(); it != MG.end(); it++){
+      k++; for (int j = 0; j < it->second.size(); j++)  ke++;
+    }
+
+    graphicg << k << " " << ke << endl; int kk = 0;
+    for (SparseGraph::iterator it = MG.begin(); it != MG.end(); it++){graphicg << kk << " " << colors[it->first] << endl; kk++;}
+    for (SparseGraph::iterator it = MG.begin(); it != MG.end(); it++)
+      for (int j = 0; j < it->second.size(); j++)
+        graphicg << it->first << " " << it->second[j] << endl;
+}
+
+void compute_bootstrapped_estimator(const vector<int> & samp, char* const & cloud, const int & Nboot, double** dist, Cloud & Cboot, \
+                                           const int & Nsubs, const double & gain){
+
+    int n = Cboot.size();
+    Cover I; MapperElements M;
+    double maxf = Cboot[0].func; double minf = Cboot[0].func;
+    for(int i = 1; i < n; i++){ maxf = max(maxf,Cboot[i].func); minf = min(minf,Cboot[i].func);}
+    int m = floor(n/pow(log(n)/log(CONSTANT),1+ETA));
+    double lip = 0; double delta = 0; vector<int> samples(m);
+    #pragma omp parallel for
+    for (int i = 0; i < Nsubs; i++){
+      SampleWithoutReplacement(n,m,samples);
+      double hausdorff_dist = 0;
+      for (int j = 0; j < n; j++){
+        Point p = Cboot[j]; double mj = dist[samp[j]][samp[samples[0]]]; double Mi = abs(p.func-Cboot[0].func)/dist[samp[j]][samp[0]];
+        for (int k = 1; k < m; k++)  mj = min(mj, dist[samp[j]][samp[samples[k]]]);
+        if(j < n-1)
+          for (int k = j+1; k < n; k++)  Mi = max(Mi, abs(p.func-Cboot[k].func)/dist[samp[j]][samp[k]]);
+        hausdorff_dist = max(hausdorff_dist, mj); lip = max(lip,Mi);
+      }
+      delta += hausdorff_dist;
+    }
+    delta /= Nsubs; double resolution;
+    if(DELTA)  resolution = lip*delta; else  resolution = lip*delta/gain;
+    //cout << delta << " " << lip << " " << resolution << endl;
+
+    AdjacencyMatrix G; G.clear(); vector<Point> adj;
+
+    for(int i = 0; i < n; i++){
+      adj.clear();
+      for(int j = 0; j < n; j++)
+        if(dist[samp[i]][samp[j]] <= delta)
+          adj.push_back(Cboot[j]);
+      G.insert(pair<Point, vector<Point> >(Cboot[i],adj));
+    }
+
+    sort(Cboot.begin(),Cboot.end());
+    I = Cover(Cboot.begin()->func, (Cboot.end()-1)->func, resolution, gain, SHIFT, 1);
+    I.proper_value();
+
+    map<int,double> colors; colors.clear(); map<int,int> num; num.clear();
+    map<Point,vector<int> > clusters; vector<int> dum; dum.clear();
+    for(int i = 0; i < Cboot.size(); i++)  clusters.insert(pair<Point,vector<int> >(Cboot[i],dum)); vector<double> col(n);
+
+    M = MapperElts(G,I,clusters,col);
+    SparseGraph MG;
+    if(DELTA)  MG = build_mapperDelta_graph(M,G,clusters,colors,num); else  MG = build_mapper_graph(M,colors,num);
+
+    char mappg[100]; sprintf(mappg,"%s_mapper_%d.txt",cloud,Nboot);
+    ofstream graphicg(mappg);
+
+    double maxv, minv; maxv = numeric_limits<double>::min(); minv = numeric_limits<double>::max();
+    for (map<int,double>::iterator iit = colors.begin(); iit != colors.end(); iit++){
+      if(iit->second > maxv)  maxv = iit->second;
+      if(minv > iit->second)  minv = iit->second;
+    }
+
+    int k = 0; int ke = 0;
+
+    for (SparseGraph::iterator it = MG.begin(); it != MG.end(); it++){
+      k++; for (int j = 0; j < it->second.size(); j++)  ke++;
+    }
+
+    graphicg << k << " " << ke << endl; int kk = 0;
+    for (SparseGraph::iterator it = MG.begin(); it != MG.end(); it++){graphicg << kk << " " << colors[it->first] << endl; kk++;}
+    for (SparseGraph::iterator it = MG.begin(); it != MG.end(); it++)
+      for (int j = 0; j < it->second.size(); j++)
+        graphicg << it->first << " " << it->second[j] << endl;
+    graphicg.close();
+
+}
+
+
+int main(int argc, char** argv){
+
+    if(argc <= 7 || argc >= 9){
+      cout << "./bottleBoot <cloud_file> <VNE> <function:/coordinate:/eccentricity:> <func_file/number/distance matrix> <gain> <Nsubs> <Nboot>" << endl;
+      return 0;}
+
+    char* const cloud = argv[1]; bool normalized = atoi(argv[2]); char* const funct = argv[3]; double gain = atof(argv[5]);
+    int Nsubs = atoi(argv[6]); int Nboot = atoi(argv[7]);
+
+    Cloud C, D, Cboot; C = read_cloud(cloud); int n = C.size();
+
+    if (strcmp(funct, "function:") == 0){
+      char* const func = argv[4]; read_function_from_file(func,C);
+    }
+    if (strcmp(funct, "coordinate:") == 0){
+      int number = atoi(argv[4]); read_coordinate(number,C);
+    }
+    if (strcmp(funct, "eccentricity") == 0){
+      char* const matrix = argv[4]; compute_eccentricity(C,matrix);
+    }
+
+    D = C;
+
+    cout << "Computing estimator..." << endl;
+    compute_estimator(C,cloud,normalized,Nsubs,gain);
+    cout << "Estimator computed." << endl;
+
+    char nam[100]; sprintf(nam,"%s_mapper.txt",cloud); char dg[100]; sprintf(dg,"%s_ExDg",cloud);
+    char command[100];
+    sprintf(command,"cat %s | ../persistence/filtgraph 1 | ../persistence/pers 1 >> %s",nam,dg); system(command);
+    char list[100]; sprintf(list,"%s_quantile",cloud);
+
+    char name[100];
+    if(normalized)  sprintf(name, "%s_VNdist", cloud);
+    else  sprintf(name, "%s_dist", cloud);
+    ifstream input(name, ios::out | ios::binary); double d;
+    double** dist = new double*[n]; for(int i = 0; i < n; i++)  dist[i] = new double[n];
+    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();
+
+    for(int i = 0; i < Nboot; i++){
+
+      Cboot.clear();
+      vector<int> samples(n); SampleWithReplacement(n,n,samples); sort(samples.begin(), samples.end());
+      for (int i = 0; i < n; i++)  Cboot.push_back(D[samples[i]]);
+
+      cout << "Computing " << i << "th bootstrapped estimator..." << endl;
+      compute_bootstrapped_estimator(samples,cloud,i,dist,Cboot,Nsubs,gain);
+      cout << "Done." << endl;
+
+      char namei[100]; sprintf(namei,"%s_mapper_%d.txt",cloud,i); char dgi[100]; sprintf(dgi,"%s_ExDg_%d",cloud,i);
+      char command[100];
+      sprintf(command,"cat %s | ../persistence/filtgraph 1 | ../persistence/pers 1 >> %s",namei,dgi); system(command);
+      sprintf(command,"../persistence/bottleneck_dist %s %s >> %s",dg,dgi,list); system(command);
+      sprintf(command, "rm %s %s",namei, dgi); system(command);
+
+    }
+
+    delete [] dist;
+
+    return 0;
+}
diff --git a/src/Nerve_GIC/example/density_smoother.cpp b/src/Nerve_GIC/example/density_smoother.cpp
new file mode 100644
index 00000000..4f3c4137
--- /dev/null
+++ b/src/Nerve_GIC/example/density_smoother.cpp
@@ -0,0 +1,155 @@
+#include <cstdlib>
+#include <string.h>
+#include <stdio.h>
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <vector>
+#include <algorithm>
+#include <set>
+#include <map>
+#include <limits>
+#include <assert.h>
+#include <math.h>
+#include <iterator>
+#include <omp.h>
+#include <list>
+#include <random>
+
+using namespace std;
+
+double EuclideanDistance(const vector<double> & V1, const vector<double> & V2){
+  assert(V1.size() == V2.size());
+  double res = 0;
+  for(int i = 0; i < V1.size(); i++)  res += pow(V1[i]-V2[i],2);
+  return sqrt(res);
+}
+
+int main(int argc, char** argv){
+
+  char* const cloud = argv[1];
+  char* const method = argv[2];
+
+  if(argc <= 4){cout << "Usage: <cloud_file> <DTM/GDE> <k/h> <threshold>" << endl; return 0;}
+  if(argc >= 6){cout << "Too many arguments !" << endl; return 0;}
+
+  ifstream Cloud(cloud);
+  vector<vector<double> > C, distances; C.clear(); double d;
+  string line;
+
+  int nb = 0; cout << "Reading cloud..." << endl;
+  while(getline(Cloud,line)){
+    if(nb%100000 == 0)  cout << "  " << nb << endl;
+    vector<double> P;
+    stringstream stream(line);
+    while(stream >> d)  P.push_back(d);
+    C.push_back(P); nb++;
+  }
+  cout << endl;
+
+  for(int i = 0; i < nb; i++){vector<double> vect(nb); distances.push_back(vect);}
+
+  char name[100]; sprintf(name,"%s_dist",cloud);
+  ifstream input(name, ios::out | ios::binary);
+
+  if(input.good()){
+
+    cout << "  Reading distances..." << endl;
+
+    for(int i = 0; i < nb; i++){
+      for (int j = 0; j < nb; j++){
+        input.read((char*) &d,8); distances[i][j] = d;
+      }
+    }
+
+    input.close();
+    cout << "  Done." << endl;
+
+  }
+
+  else{
+
+    cout << "  Computing distances..." << endl;
+    input.close(); ofstream output(name, ios::out | ios::binary);
+
+    for(int i = 0; i < nb-1; i++){
+      if( (int) floor( 100*((double) i)/((double) nb)+1 ) %10 == 0  ){cout << "\r" << floor( 100*((double) i)/((double) nb) +1) << "%" << flush;}
+      for (int j = i+1; j < nb; j++){
+        d = EuclideanDistance(C[i],C[j]); distances[i][j] = d; distances[j][i] = d;
+        output.write((char*) &d,8);
+      }
+    }
+
+    output.close();
+    cout << endl << "  Done." << endl;
+
+  }
+
+  if(strcmp(method,"DTM") == 0){
+
+    int k = atoi(argv[3]);
+    double thresh = atof(argv[4]);
+
+    char nameo[100]; sprintf(nameo,"%s_DTM_%d_%g",cloud,k,thresh);
+    ofstream output(nameo); vector<double> densities; double M = 0;
+
+    cout << "Computing DTM..." << endl;
+    for(int i = 0; i < nb; i++){
+      if( (int) floor( 100*((double) i)/((double) nb)+1 ) %10 == 0  ){cout << "\r" << floor( 100*((double) i)/((double) nb) +1) << "%" << flush;}
+      double density = 0;
+      sort(distances[i].begin(),distances[i].end());
+      for(int j = 0; j < k; j++)
+        density += pow(distances[i][j],2);
+      density = sqrt(k)/sqrt(density);
+      if(density >= M)  M = density;
+      densities.push_back(density);
+
+    }
+    cout << endl << "  Done." << endl;
+
+    for(int i = 0; i < nb; i++){
+      if(densities[i] >= thresh*M){
+        for(int j = 0; j < C[i].size(); j++)  output << C[i][j] << " ";
+        output << endl;
+      }
+    }
+
+    output.close();
+
+  }
+
+  if(strcmp(method,"GaussianDE") == 0){
+
+    double bandwidth = atof(argv[3]);
+    double thresh = atof(argv[4]);
+
+    char nameo[100]; sprintf(nameo,"%s_GDE_%g_%g",cloud,bandwidth,thresh);
+    ofstream output(nameo); vector<double> densities; double M = 0;
+
+    cout << "Computing GDE..." << endl;
+    for(int i = 0; i < nb; i++){
+      if( (int) floor( 100*((double) i)/((double) nb)+1 ) %10 == 0  ){cout << "\r" << floor( 100*((double) i)/((double) nb) +1) << "%" << flush;}
+      double density = 0;
+      for(int j = 0; j < nb; j++)
+        density += exp(-pow(distances[i][j],2)/(2*pow(bandwidth,2)));
+      density /= sqrt(2*3.1415)*nb*bandwidth;
+      if(density >= M)  M = density;
+      densities.push_back(density);
+
+    }
+    cout << endl << "  Done." << endl;
+
+    for(int i = 0; i < nb; i++){
+      if(densities[i] >= thresh*M){
+        for(int j = 0; j < C[i].size(); j++)  output << C[i][j] << " ";
+        output << endl;
+      }
+    }
+
+    output.close();
+
+  }
+
+  return 0;
+
+}
diff --git a/src/Nerve_GIC/example/graph_off.cpp b/src/Nerve_GIC/example/graph_off.cpp
new file mode 100644
index 00000000..574d9fc6
--- /dev/null
+++ b/src/Nerve_GIC/example/graph_off.cpp
@@ -0,0 +1,49 @@
+#include <iostream>
+#include <set>
+#include <fstream>
+#include <vector>
+#include <string.h>
+#include <sstream>
+#include <stdlib.h>
+#include <algorithm>
+
+using namespace std;
+
+int main (int argc, char **argv) {
+
+  char* const fileoff = argv[1]; ifstream input(fileoff); char ngoff[100]; sprintf(ngoff,"%s_NG", fileoff); ofstream output(ngoff);
+  char buf[256];
+  vector<vector<int> > NG;
+
+  input.getline(buf, 255);  // skip "OFF"
+  int n, m;
+  input >> n >> m;
+  input.getline(buf, 255);  // skip "0"
+
+  // read vertices
+  double x,y,z; vector<int> ng; int nn = n;
+  while(nn-->0) {
+    input >> x >> z >> y; NG.push_back(ng);
+  }
+
+  // read triangles
+  int d, p, q, s;
+  while (m-->0) {
+    input >> d >> p >> q >> s;
+    NG[p].push_back(q); NG[p].push_back(s);
+    NG[q].push_back(p); NG[q].push_back(s);
+    NG[s].push_back(q); NG[s].push_back(p);
+  }
+
+  for(int i = 0; i < n; i++){
+    vector<int> ng = NG[i];
+    sort(ng.begin(),ng.end()); vector<int>::iterator iter = unique(ng.begin(),ng.end()); ng.resize(distance(ng.begin(),iter));
+    int size = ng.size();
+    for(int j = 0; j < size; j++)
+      output << ng[j] << " ";
+    output << endl;
+  }
+
+  return 0;
+
+}
diff --git a/src/Nerve_GIC/example/km.py b/src/Nerve_GIC/example/km.py
new file mode 100755
index 00000000..d3cfad59
--- /dev/null
+++ b/src/Nerve_GIC/example/km.py
@@ -0,0 +1,390 @@
+from __future__ import division
+import numpy as np
+from collections import defaultdict
+import json
+import itertools
+from sklearn import cluster, preprocessing, manifold
+from datetime import datetime
+import sys
+
+class KeplerMapper(object):
+  # With this class you can build topological networks from (high-dimensional) data.
+  #
+  # 1)   	Fit a projection/lens/function to a dataset and transform it. 
+  #     	For instance "mean_of_row(x) for x in X"
+  # 2)   	Map this projection with overlapping intervals/hypercubes. 
+  #    		Cluster the points inside the interval 
+  #    		(Note: we cluster on the inverse image/original data to lessen projection loss).
+  #    		If two clusters/nodes have the same members (due to the overlap), then: 
+  #    		connect these with an edge.
+  # 3)  	Visualize the network using HTML and D3.js.
+  # 
+  # functions
+  # ---------
+  # fit_transform:   Create a projection (lens) from a dataset
+  # map:         	Apply Mapper algorithm on this projection and build a simplicial complex
+  # visualize:    	Turns the complex dictionary into a HTML/D3.js visualization
+  
+  def __init__(self, verbose=2):
+    self.verbose = verbose
+    
+    self.chunk_dist = []
+    self.overlap_dist = []
+    self.d = []
+    self.nr_cubes = 0
+    self.overlap_perc = 0
+    self.clusterer = False
+
+  def fit_transform(self, X, projection="sum", scaler=preprocessing.MinMaxScaler()):
+    # Creates the projection/lens from X. 
+    #
+    # Input:      X. Input features as a numpy array.
+    # Output:     projected_X. original data transformed to a projection (lens).
+    # 
+    # parameters
+    # ----------
+    # projection:   Projection parameter is either a string, 
+    #               a scikit class with fit_transform, like manifold.TSNE(), 
+    #               or a list of dimension indices.
+    # scaler:       if None, do no scaling, else apply scaling to the projection
+    #               Default: Min-Max scaling
+    
+    self.scaler = scaler
+    self.projection = str(projection)
+    
+    # Detect if projection is a class (for scikit-learn)
+    if str(type(projection))[1:6] == "class": #TODO: de-ugly-fy
+      reducer = projection
+      if self.verbose > 0:
+        try:    
+          projection.set_params(**{"verbose":self.verbose})
+        except:
+          pass
+        print("\n..Projecting data using: \n\t%s\n"%str(projection))
+      X = reducer.fit_transform(X)
+    
+    # Detect if projection is a string (for standard functions)
+    if isinstance(projection, str):
+      if self.verbose > 0:
+        print("\n..Projecting data using: %s"%(projection))
+      # Stats lenses
+      if projection == "sum": # sum of row
+        X = np.sum(X, axis=1).reshape((X.shape[0],1))
+      if projection == "mean": # mean of row
+        X = np.mean(X, axis=1).reshape((X.shape[0],1))
+      if projection == "median": # mean of row
+        X = np.median(X, axis=1).reshape((X.shape[0],1))
+      if projection == "max": # max of row
+        X = np.max(X, axis=1).reshape((X.shape[0],1))
+      if projection == "min": # min of row
+        X = np.min(X, axis=1).reshape((X.shape[0],1))
+      if projection == "std": # std of row
+        X = np.std(X, axis=1).reshape((X.shape[0],1))
+        
+      if projection == "dist_mean": # Distance of x to mean of X
+        X_mean = np.mean(X, axis=0) 
+        X = np.sum(np.sqrt((X - X_mean)**2), axis=1).reshape((X.shape[0],1))
+
+    # Detect if projection is a list (with dimension indices)
+    if isinstance(projection, list):
+      if self.verbose > 0:
+        print("\n..Projecting data using: %s"%(str(projection)))
+      X = X[:,np.array(projection)]
+      
+    # Scaling
+    if scaler is not None:
+      if self.verbose > 0:
+        print("\n..Scaling with: %s\n"%str(scaler))
+      X = scaler.fit_transform(X)
+    
+    return X
+
+  def map(self, projected_X, inverse_X=None, clusterer=cluster.DBSCAN(eps=0.5,min_samples=3), nr_cubes=10, overlap_perc=0.1):
+    # This maps the data to a simplicial complex. Returns a dictionary with nodes and links.
+    # 
+    # Input:    projected_X. A Numpy array with the projection/lens. 
+    # Output:    complex. A dictionary with "nodes", "links" and "meta information"
+    #
+    # parameters
+    # ----------
+    # projected_X  	projected_X. A Numpy array with the projection/lens. Required.
+    # inverse_X    	Numpy array or None. If None then the projection itself is used for clustering.
+    # clusterer    	Scikit-learn API compatible clustering algorithm. Default: DBSCAN
+    # nr_cubes    	Int. The number of intervals/hypercubes to create.
+    # overlap_perc  Float. The percentage of overlap "between" the intervals/hypercubes.
+    
+    start = datetime.now()
+    
+    # Helper function
+    def cube_coordinates_all(nr_cubes, nr_dimensions):
+      # Helper function to get origin coordinates for our intervals/hypercubes
+      # Useful for looping no matter the number of cubes or dimensions
+      # Example:   	if there are 4 cubes per dimension and 3 dimensions 
+      #       		return the bottom left (origin) coordinates of 64 hypercubes, 
+      #       		as a sorted list of Numpy arrays
+      # TODO: elegance-ify...
+      l = []
+      for x in range(nr_cubes):
+        l += [x] * nr_dimensions
+      return [np.array(list(f)) for f in sorted(set(itertools.permutations(l,nr_dimensions)))]
+    
+    nodes = defaultdict(list)
+    links = defaultdict(list)
+    complex = {}
+    self.nr_cubes = nr_cubes
+    self.clusterer = clusterer
+    self.overlap_perc = overlap_perc
+    
+    if self.verbose > 0:
+      print("Mapping on data shaped %s using dimensions\n"%(str(projected_X.shape)))
+    
+    # If inverse image is not provided, we use the projection as the inverse image (suffer projection loss)
+    if inverse_X is None:
+      inverse_X = projected_X
+      
+    # We chop up the min-max column ranges into 'nr_cubes' parts
+    self.chunk_dist = (np.max(projected_X, axis=0) - np.min(projected_X, axis=0))/nr_cubes
+
+    # We calculate the overlapping windows distance 
+    self.overlap_dist = self.overlap_perc * self.chunk_dist
+
+    # We find our starting point
+    self.d = np.min(projected_X, axis=0)
+    
+    # Use a dimension index array on the projected X
+    # (For now this uses the entire dimensionality, but we keep for experimentation)
+    di = np.array([x for x in range(projected_X.shape[1])])
+    
+    # Prefix'ing the data with ID's
+    ids = np.array([x for x in range(projected_X.shape[0])])
+    projected_X = np.c_[ids,projected_X]
+    inverse_X = np.c_[ids,inverse_X]
+
+    # Subdivide the projected data X in intervals/hypercubes with overlap
+    if self.verbose > 0:
+      total_cubes = len(cube_coordinates_all(nr_cubes,projected_X.shape[1]))
+      print("Creating %s hypercubes."%total_cubes)
+
+    for i, coor in enumerate(cube_coordinates_all(nr_cubes,di.shape[0])):
+      # Slice the hypercube
+      hypercube = projected_X[ np.invert(np.any((projected_X[:,di+1] >= self.d[di] + (coor * self.chunk_dist[di])) & 
+          (projected_X[:,di+1] < self.d[di] + (coor * self.chunk_dist[di]) + self.chunk_dist[di] + self.overlap_dist[di]) == False, axis=1 )) ]
+      
+      if self.verbose > 1:
+        print("There are %s points in cube_%s / %s with starting range %s"%
+              (hypercube.shape[0],i,total_cubes,self.d[di] + (coor * self.chunk_dist[di])))
+      
+      # If at least one sample inside the hypercube
+      if hypercube.shape[0] > 0:
+        # Cluster the data point(s) in the cube, skipping the id-column
+        # Note that we apply clustering on the inverse image (original data samples) that fall inside the cube.
+        inverse_x = inverse_X[[int(nn) for nn in hypercube[:,0]]]
+        
+        clusterer.fit(inverse_x[:,1:])
+        
+        if self.verbose > 1:
+          print("Found %s clusters in cube_%s\n"%(np.unique(clusterer.labels_[clusterer.labels_ > -1]).shape[0],i))
+        
+        #Now for every (sample id in cube, predicted cluster label)
+        for a in np.c_[hypercube[:,0],clusterer.labels_]:
+          if a[1] != -1: #if not predicted as noise
+            cluster_id = str(coor[0])+"_"+str(i)+"_"+str(a[1])+"_"+str(coor)+"_"+str(self.d[di] + (coor * self.chunk_dist[di])) # TODO: de-rudimentary-ify
+            nodes[cluster_id].append( int(a[0]) ) # Append the member id's as integers
+      else:
+        if self.verbose > 1:
+          print("Cube_%s is empty.\n"%(i))
+
+    # Create links when clusters from different hypercubes have members with the same sample id.
+    candidates = itertools.combinations(nodes.keys(),2)
+    for candidate in candidates:
+      # if there are non-unique members in the union
+      if len(nodes[candidate[0]]+nodes[candidate[1]]) != len(set(nodes[candidate[0]]+nodes[candidate[1]])):
+        links[candidate[0]].append( candidate[1] )
+
+    # Reporting
+    if self.verbose > 0:
+      nr_links = 0
+      for k in links:
+        nr_links += len(links[k])
+      print("\ncreated %s edges and %s nodes in %s."%(nr_links,len(nodes),str(datetime.now()-start)))
+    
+    complex["nodes"] = nodes
+    complex["links"] = links
+    complex["meta"] = self.projection
+
+    return complex
+
+  def visualize(self, complex, color_function="", path_html="mapper_visualization_output.html", title="My Data", 
+          graph_link_distance=30, graph_gravity=0.1, graph_charge=-120, custom_tooltips=None, width_html=0, 
+          height_html=0, show_tooltips=True, show_title=True, show_meta=True, res=0,gain=0,minimum=0,maximum=0):
+    # Turns the dictionary 'complex' in a html file with d3.js
+    #
+    # Input:      complex. Dictionary (output from calling .map())
+    # Output:      a HTML page saved as a file in 'path_html'.
+    # 
+    # parameters
+    # ----------
+    # color_function    	string. Not fully implemented. Default: "" (distance to origin)
+    # path_html        		file path as string. Where to save the HTML page.
+    # title          		string. HTML page document title and first heading.
+    # graph_link_distance  	int. Edge length.
+    # graph_gravity     	float. "Gravity" to center of layout.
+    # graph_charge      	int. charge between nodes.
+    # custom_tooltips   	None or Numpy Array. You could use "y"-label array for this.
+    # width_html        	int. Width of canvas. Default: 0 (full width)
+    # height_html       	int. Height of canvas. Default: 0 (full height)
+    # show_tooltips     	bool. default:True
+    # show_title      		bool. default:True
+    # show_meta        		bool. default:True
+
+    # Format JSON for D3 graph
+    json_s = {}
+    json_s["nodes"] = []
+    json_s["links"] = []
+    k2e = {} # a key to incremental int dict, used for id's when linking
+
+    for e, k in enumerate(complex["nodes"]):
+      # Tooltip and node color formatting, TODO: de-mess-ify
+      if custom_tooltips is not None:
+        tooltip_s = "<h2>Cluster %s</h2>"%k + " ".join(str(custom_tooltips[complex["nodes"][k][0]]).split(" "))
+        if color_function == "Lens average":
+          tooltip_i = int(30*(custom_tooltips[complex["nodes"][k][0]]-minimum)/(maximum-minimum))
+          json_s["nodes"].append({"name": str(k), "tooltip": tooltip_s, "group": 2 * int(np.log(complex["nodes"][k][2])), "color": tooltip_i})
+        else:
+          json_s["nodes"].append({"name": str(k), "tooltip": tooltip_s, "group": 2 * int(np.log(len(complex["nodes"][k]))), "color": str(k.split("_")[0])})
+      else:
+        tooltip_s = "<h2>Cluster %s</h2>Contains %s members."%(k,len(complex["nodes"][k]))
+        json_s["nodes"].append({"name": str(k), "tooltip": tooltip_s, "group": 2 * int(np.log(len(complex["nodes"][k]))), "color": str(k.split("_")[0])})
+      k2e[k] = e
+    for k in complex["links"]:
+      for link in complex["links"][k]:
+        json_s["links"].append({"source": k2e[k], "target":k2e[link],"value":1})
+
+    # Width and height of graph in HTML output
+    if width_html == 0:
+      width_css = "100%"
+      width_js = 'document.getElementById("holder").offsetWidth-20'
+    else:
+      width_css = "%spx" % width_html
+      width_js = "%s" % width_html
+    if height_html == 0:
+      height_css = "100%"
+      height_js = 'document.getElementById("holder").offsetHeight-20'
+    else:
+      height_css = "%spx" % height_html
+      height_js = "%s" % height_html
+    
+    # Whether to show certain UI elements or not
+    if show_tooltips == False:
+      tooltips_display = "display: none;"
+    else:
+      tooltips_display = ""
+      
+    if show_meta == False:
+      meta_display = "display: none;"
+    else:
+      meta_display = ""
+      
+    if show_title == False:
+      title_display = "display: none;"
+    else:
+      title_display = ""  
+    
+    with open(path_html,"wb") as outfile:
+      html = """<!DOCTYPE html>
+    <meta charset="utf-8">
+    <meta name="generator" content="KeplerMapper">
+    <title>%s | KeplerMapper</title>
+    <link href='https://fonts.googleapis.com/css?family=Roboto:700,300' rel='stylesheet' type='text/css'>
+    <style>
+    * {margin: 0; padding: 0;}
+    html { height: 100%%;}
+    body {background: #111; height: 100%%; font: 100 16px Roboto, Sans-serif;}
+    .link { stroke: #999; stroke-opacity: .333;  }
+    .divs div { border-radius: 50%%; background: red; position: absolute; }
+    .divs { position: absolute; top: 0; left: 0; }
+    #holder { position: relative; width: %s; height: %s; background: #111; display: block;}
+    h1 { %s padding: 20px; color: #fafafa; text-shadow: 0px 1px #000,0px -1px #000; position: absolute; font: 300 30px Roboto, Sans-serif;}
+    h2 { text-shadow: 0px 1px #000,0px -1px #000; font: 700 16px Roboto, Sans-serif;}
+    .meta {  position: absolute; opacity: 0.9; width: 220px; top: 80px; left: 20px; display: block; %s background: #000; line-height: 25px; color: #fafafa; border: 20px solid #000; font: 100 16px Roboto, Sans-serif;}
+    div.tooltip { position: absolute; width: 380px; display: block; %s padding: 20px; background: #000; border: 0px; border-radius: 3px; pointer-events: none; z-index: 999; color: #FAFAFA;}
+    }
+    </style>
+    <body>
+    <div id="holder">
+      <h1>%s</h1>
+      <p class="meta">
+      <b>Lens</b><br>%s<br><br>
+      <b>Length of intervals</b><br>%s<br><br>
+      <b>Overlap percentage</b><br>%s%%<br><br>
+      <b>Color Function</b><br>%s
+      </p>
+    </div>
+    <script src="https://cdnjs.cloudflare.com/ajax/libs/d3/3.5.5/d3.min.js"></script>
+    <script>
+    var width = %s,
+      height = %s;
+    var color = d3.scale.ordinal()
+      .domain(["0","1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13","14","15","16","17","18","19","20","21","22","23","24","25","26","27","28","29","30"])
+      .range(["#FF0000","#FF1400","#FF2800","#FF3c00","#FF5000","#FF6400","#FF7800","#FF8c00","#FFa000","#FFb400","#FFc800","#FFdc00","#FFf000","#fdff00","#b0ff00","#65ff00","#17ff00","#00ff36","#00ff83","#00ffd0","#00e4ff","#00c4ff","#00a4ff","#00a4ff","#0084ff","#0064ff","#0044ff","#0022ff","#0002ff","#0100ff","#0300ff","#0500ff"]);
+    var force = d3.layout.force()
+      .charge(%s)
+      .linkDistance(%s)
+      .gravity(%s)
+      .size([width, height]);
+    var svg = d3.select("#holder").append("svg")
+      .attr("width", width)
+      .attr("height", height);
+    
+    var div = d3.select("#holder").append("div")   
+      .attr("class", "tooltip")               
+      .style("opacity", 0.0);
+    
+    var divs = d3.select('#holder').append('div')
+      .attr('class', 'divs')
+      .attr('style', function(d) { return 'overflow: hidden; width: ' + width + 'px; height: ' + height + 'px;'; });  
+    
+      graph = %s;
+      force
+        .nodes(graph.nodes)
+        .links(graph.links)
+        .start();
+      var link = svg.selectAll(".link")
+        .data(graph.links)
+        .enter().append("line")
+        .attr("class", "link")
+        .style("stroke-width", function(d) { return Math.sqrt(d.value); });
+      var node = divs.selectAll('div')
+      .data(graph.nodes)
+        .enter().append('div')
+        .on("mouseover", function(d) {      
+          div.transition()        
+            .duration(200)      
+            .style("opacity", .9);
+          div .html(d.tooltip + "<br/>")  
+            .style("left", (d3.event.pageX + 100) + "px")     
+            .style("top", (d3.event.pageY - 28) + "px");    
+          })                  
+        .on("mouseout", function(d) {       
+          div.transition()        
+            .duration(500)      
+            .style("opacity", 0);   
+        })
+        .call(force.drag);
+      
+      node.append("title")
+        .text(function(d) { return d.name; });
+      force.on("tick", function() {
+      link.attr("x1", function(d) { return d.source.x; })
+        .attr("y1", function(d) { return d.source.y; })
+        .attr("x2", function(d) { return d.target.x; })
+        .attr("y2", function(d) { return d.target.y; });
+      node.attr("cx", function(d) { return d.x; })
+        .attr("cy", function(d) { return d.y; })
+        .attr('style', function(d) { return 'width: ' + (d.group * 2) + 'px; height: ' + (d.group * 2) + 'px; ' + 'left: '+(d.x-(d.group))+'px; ' + 'top: '+(d.y-(d.group))+'px; background: '+color(d.color)+'; box-shadow: 0px 0px 3px #111; box-shadow: 0px 0px 33px '+color(d.color)+', inset 0px 0px 5px rgba(0, 0, 0, 0.2);'})
+        ;
+      });
+    </script>"""%(title,width_css, height_css, title_display, meta_display, tooltips_display, title,complex["meta"],res,gain*100,color_function,width_js,height_js,graph_charge,graph_link_distance,graph_gravity,json.dumps(json_s))
+      outfile.write(html.encode("utf-8"))
+    if self.verbose > 0:
+      print("\nWrote d3.js graph to '%s'"%path_html)
\ No newline at end of file
diff --git a/src/Nerve_GIC/example/mapper.cpp b/src/Nerve_GIC/example/mapper.cpp
new file mode 100755
index 00000000..05dd87dd
--- /dev/null
+++ b/src/Nerve_GIC/example/mapper.cpp
@@ -0,0 +1,179 @@
+#include "mapper.h"
+#define SHIFT 1
+#define DELTA 1
+#define RESOLUTION 1
+
+int main(int argc, char** argv){
+
+  if(argc <= 10 || argc >= 12){cout << "./mapper <cloud_file> <function:/coordinate:/eccentricity:> <func_file/number/matrix_file> <graph:/parameter:/percentage:>" <<\
+                                       " <graph_file/delta> <VNE> <cover:/uniform:> <cover_file/resolution> <gain> <mask>" << endl; return 0;}
+
+  char* const cloud = argv[1];
+  char* const funct = argv[2];
+  bool normalized = atoi(argv[6]);
+  char* const graph = argv[4];
+  int mask;
+  char* const covering = argv[7];
+
+  Cover I; AdjacencyMatrix G; Cloud C;
+  MapperElements M;
+
+  cout << "Reading input cloud from file " << cloud << "..." << endl;
+  C = read_cloud(cloud);
+  cout << "  Done." << endl;
+
+  double r,g; char namefunc[100];
+
+  if (strcmp(funct, "function:") == 0){
+    char* const func = argv[3];
+    cout << "Reading input filter from file " << func << "..." << endl;
+    read_function_from_file(func,C); sprintf(namefunc,"%s",func);
+  }
+  if (strcmp(funct, "coordinate:") == 0){
+    int number = atoi(argv[3]);
+    cout << "Using coordinate " << number << " as a filter..." << endl;
+    read_coordinate(number,C); sprintf(namefunc,"Coordinate %d",number);
+  }
+  if (strcmp(funct, "eccentricity:") == 0){
+    char* const matrix = argv[3];
+    cout << "Computing eccentricity with distance matrix " << matrix << "..." << endl;
+    compute_eccentricity(C,matrix); sprintf(namefunc,"eccentricity");
+  }
+  cout << "  Done." << endl;
+
+  if (strcmp(graph, "graph:") == 0){
+    char* const graph_name = argv[5];
+    cout << "Reading neighborhood graph from file " << graph_name << "..." << endl;
+    G = build_neighborhood_graph_from_file(C,graph_name);
+  }
+  if (strcmp(graph, "percentage:") == 0){
+    double delta = atof(argv[5]);
+    char name1[100]; sprintf(name1, "%s_dist", cloud); char* const name2 = name1;
+    cout << "Computing neighborhood graph with delta percentage = " << delta << "..." << endl;
+    G = build_neighborhood_graph_from_scratch_with_percentage(C,delta,name2,normalized);
+  }
+  if (strcmp(graph, "parameter:") == 0){
+    double delta = atof(argv[5]);
+    char name1[100]; sprintf(name1, "%s_dist", cloud); char* const name2 = name1;
+    cout << "Computing neighborhood graph with delta parameter = " << delta << "..." << endl;
+    G = build_neighborhood_graph_from_scratch_with_parameter(C,delta,name2,normalized);
+  }
+  cout << "  Done." << endl;
+
+  cout << "Sorting cloud..." << endl;
+  sort(C.begin(),C.end());
+  cout << "  Done." << endl;
+
+  if(strcmp(covering,"cover:") == 0){
+    char* const cover = argv[8]; mask = atoi(argv[9]);
+    cout << "Reading user-defined cover from file " << cover << "..." << endl;
+    I = Cover(cover); I.sort_covering();
+    assert (I.intervals[0].first <= C.begin()->func && I.intervals[I.intervals.size()-1].second >= (C.end()-1)->func);
+  }
+  else{
+    double resolution = atof(argv[8]); r = resolution;
+    double gain = atof(argv[9]); mask = atoi(argv[10]); g = gain;
+    cout << "Computing uniform cover with resolution " << resolution << " and gain " << gain << "..." << endl;
+    I = Cover(C.begin()->func, (C.end()-1)->func, resolution, gain, SHIFT, RESOLUTION);
+  }
+  I.proper_value();
+  /*for (int i = 0; i < I.res; i++)
+    cout << "  " << I.intervals[i].first << " " << I.intervals[i].second << " " << I.value[i] << endl;
+  cout << "  Done." << endl;*/
+
+  map<Point,vector<int> > clusters; vector<int> dum; dum.clear();
+  for(int i = 0; i < G.size(); i++)  clusters.insert(pair<Point,vector<int> >(C[i],dum)); vector<double> col(G.size());
+  cout << "Computing Mapper nodes..." << endl;
+  M = MapperElts(G,I,clusters,col);
+  cout << "  Done." << endl;
+
+  cout << "Computing intersections..." << endl;
+  map<int,double> colors; colors.clear(); map<int,int> num; num.clear(); SparseGraph MG;
+  if(!DELTA)  MG = build_mapper_graph(M,colors,num);
+  else  MG = build_mapperDelta_graph(M,G,clusters,colors,num);
+  cout << "  Done." << endl;
+
+  /*cout << "Computing Nerve..." << endl;
+  SimplicialComplex Nerve = compute_Nerve(C,clusters);
+  for(int i = 0; i < Nerve.size(); i++){
+    for(int j = 0; j < Nerve[i].size(); j++)  cout << Nerve[i][j] << " ";
+    cout << endl;
+  }*/
+
+  /*cout << "Computing GIC..." << endl;
+  SimplicialComplex GIC = compute_GIC(G,clusters);
+  for(int i = 0; i < GIC.size(); i++){
+    for(int j = 0; j < GIC[i].size(); j++)  cout << GIC[i][j] << " ";
+    cout << endl;
+  }*/
+
+  //cout << "Checking intersection crossings..." << endl;
+  //vector<pair<int,int> > error_pairs = check_intersection_crossing(M,G);
+  //cout << "Done." << endl;
+
+  cout << "Writing outputs..." << endl;
+
+  char mapp[11] = "mapper.dot";
+  char mappg[11] = "mapper.txt";
+  ofstream graphic(mapp); graphic << "graph Mapper {" << endl;
+  ofstream graphicg(mappg);
+
+  double maxv, minv; maxv = numeric_limits<double>::min(); minv = numeric_limits<double>::max();
+  for (map<int,double>::iterator iit = colors.begin(); iit != colors.end(); iit++){
+    if(iit->second > maxv){maxv = iit->second;}
+    if(minv > iit->second){minv = iit->second;}
+  }
+
+  int k = 0; vector<int> nodes; nodes.clear();
+  for (SparseGraph::iterator iit = MG.begin(); iit != MG.end(); iit++){
+    if(num[iit->first] > mask){
+      nodes.push_back(iit->first);
+      graphic << iit->first << "[shape=circle fontcolor=black color=black label=\""
+              << iit->first /*<< ":" << num[iit->first]*/ << "\" style=filled fillcolor=\""
+              << (1-(maxv-colors[iit->first])/(maxv-minv))*0.6 << ", 1, 1\"]" << endl;
+      k++;
+    }
+  }
+  int ke = 0;
+  for (SparseGraph::iterator it = MG.begin(); it != MG.end(); it++)
+    for (int j = 0; j < it->second.size(); j++)
+      if(num[it->first] > mask && num[it->second[j]] > mask){
+        graphic << "  " << it->first << " -- " << it->second[j] << " [weight=15];" << endl; ke++;}
+  graphic << "}"; graphic.close();
+
+  graphicg << cloud << endl;
+  graphicg << namefunc << endl;
+  graphicg << r << " " << g << endl;
+  graphicg << k << " " << ke << endl; int kk = 0;
+  for (vector<int>::iterator iit = nodes.begin(); iit != nodes.end(); iit++){
+    graphicg << kk << " " << colors[*iit] << " " << num[*iit] << endl; kk++;}
+  for (SparseGraph::iterator it = MG.begin(); it != MG.end(); it++)
+    for (int j = 0; j < it->second.size(); j++)
+      if(num[it->first] > mask && num[it->second[j]] > mask)
+        graphicg << it->first << " " << it->second[j] << endl;
+  graphicg.close();
+
+  cout << "  Done." << endl;
+
+  char command[100];
+  sprintf(command,"neato %s -Tpdf -o mapper.pdf",mapp); system(command);
+  sprintf(command,"python visu.py"); system(command);
+  sprintf(command,"firefox mapper_visualization_output.html"); system(command);
+  sprintf(command,"evince mapper.pdf"); system(command);
+  sprintf(command,"rm %s %s mapper_visualization_output.html mapper.pdf",mapp,mappg); system(command);
+
+  /*
+  if (error_pairs.size() > 0)
+    for (vector<pair<int,int> >::iterator it = error_pairs.begin(); it != error_pairs.end(); it++)
+      if( find(nodes.begin(),nodes.end(),it->first)!=nodes.end() && find(nodes.begin(),nodes.end(),it->second)!=nodes.end() )
+        graphic << "  " << it->first << " -- " << it->second << " [weight=15];" << endl;
+  */
+
+  int dim = C[0].coord.size();
+  if (dim <= 3){
+    plotNG(dim,G);
+    plotCover(dim,G,col,I);
+  }
+
+return 0;
+}
diff --git a/src/Nerve_GIC/example/mapper_parameter.cpp b/src/Nerve_GIC/example/mapper_parameter.cpp
new file mode 100644
index 00000000..2302d164
--- /dev/null
+++ b/src/Nerve_GIC/example/mapper_parameter.cpp
@@ -0,0 +1,165 @@
+#include "mapper.h"
+#define ETA 0.001
+#define VERBOSE 1
+#define DELTA 1
+#define RESOLUTION 1
+#define CONSTANT 10
+
+double GetUniform(){
+  static default_random_engine re;
+  static uniform_real_distribution<double> Dist(0,1);
+  return Dist(re);
+}
+
+void SampleWithoutReplacement(int populationSize, int sampleSize, vector<int> & samples){
+  int& n = sampleSize; int& N = populationSize;
+  int t = 0; int m = 0; double u;
+  while (m < n){
+    u = GetUniform();
+    if ( (N - t)*u >= n - m )
+      t++;
+    else{samples[m] = t; t++; m++;}
+  }
+}
+
+int main(int argc, char** argv){
+
+    if(argc <= 7 || argc >= 9){cout << "./param <cloud_file> <function:/coordinate:/eccentricity:> <func_file/number/matrix_file>" << \
+                                        " <VNE> <gain> <subsampling:/graph:> <N/graph_name>" << endl; return 0;}
+
+    char* const cloud = argv[1];
+    char* const funct = argv[2];
+    bool VNE = atoi(argv[4]);
+    double g = atof(argv[5]);
+    char* const method = argv[6];
+
+    Cloud C;
+
+    if(VERBOSE)  cout << "Reading input cloud from file " << cloud << "..." << endl;
+    C = read_cloud(cloud); int n = C.size();
+    if(VERBOSE)  cout << "  Done." << endl;
+
+    if (strcmp(funct, "function:") == 0){
+      char* const func = argv[3];
+      if(VERBOSE)  cout << "Reading input filter from file " << func << "..." << endl;
+      read_function_from_file(func,C);
+    }
+    if (strcmp(funct, "coordinate:") == 0){
+      int number = atoi(argv[3]);
+      if(VERBOSE)  cout << "Using coordinate " << number << " as a filter..." << endl;
+      read_coordinate(number,C);
+    }
+    if (strcmp(funct, "eccentricity:") == 0){
+      char* const matrix = argv[3];
+      cout << "Computing eccentricity with distance matrix " << matrix << "..." << endl;
+      compute_eccentricity(C,matrix);
+    }
+    if(VERBOSE)  cout << "  Done." << endl;
+
+    double maxf = C[0].func; double minf = C[0].func;
+    for(int i = 1; i < n; i++){ maxf = max(maxf,C[i].func); minf = min(minf,C[i].func);}
+
+    double delta = 0; double lip = 0;
+
+    if (strcmp(method, "subsampling:") == 0){
+
+      char name[100];
+      if(VNE)  sprintf(name,"%s_VNdist",cloud);
+      else  sprintf(name,"%s_dist",cloud);
+      ifstream input(name, ios::out | ios::binary);
+      vector<vector<double> > dist; dist.clear(); double d;
+
+      if(input.good()){
+        cout << "Reading distances..." << endl;
+        for(int i = 0; i < n; i++){
+          vector<double> dis; dis.clear();
+          for (int j = 0; j < n; j++){input.read((char*) &d,8); dis.push_back(d);}
+          dist.push_back(dis);
+        }
+        input.close();
+      }
+      else{
+        cout << "Computing distances..." << endl; vector<double> V;
+        input.close(); ofstream output(name, ios::out | ios::binary);
+        if(VNE){
+          int dim = C[0].coord.size();
+          for(int i = 0; i < dim; i++){
+            double meani = 0;
+            for (int j = 0; j < n; j++)  meani += C[j].coord[i]/n;
+            double vari = 0;
+            for (int j = 0; j < n; j++)  vari += pow(C[j].coord[i]-meani,2)/n;
+            V.push_back(vari);
+          }
+        }
+
+        for(int i = 0; i < n; i++){
+          if( (int) floor( 100*((double) i)/((double) n)+1 ) %10 == 0  ){cout << "\r" << floor( 100*((double) i)/((double) n) +1) << "%" << flush;}
+          vector<double> dis; dis.clear();
+          for (int j = 0; j < n; j++){
+            if(!VNE){d = C[i].EuclideanDistance(C[j]); dis.push_back(d);}
+            else{d = C[i].VarianceNormalizedEuclideanDistance(C[j],V); dis.push_back(d);}
+            output.write((char*) &d,8);
+          }
+          dist.push_back(dis);
+        }
+        output.close();
+      }
+
+      cout << "Done." << endl;
+
+      int m = floor(n/pow(log(n)/log(CONSTANT),1+ETA)); m = min(m,n-1);
+
+      if(VERBOSE)  cout << "dimension = " << C[0].coord.size() << endl << "n = " << n << " and s(n) = " << m << endl;
+      if(VERBOSE)  cout << "range = [" << minf << ", " << maxf << "]" << endl;
+
+      vector<int> samples(m); int N = atoi(argv[7]);
+      #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++){
+          Point p = C[j];
+          double mj = dist[j][samples[0]];
+          double Mi = abs(p.func-C[0].func)/(dist[j][0]);
+          for (int k = 1; k < m; k++){mj = min(mj, dist[j][samples[k]]);}
+          for (int k = j+1; k < n; k++){Mi = max(Mi, abs(p.func-C[k].func)/(dist[j][k]));}
+          hausdorff_dist = max(hausdorff_dist, mj); lip = max(lip,Mi);
+        }
+        delta += hausdorff_dist;
+      }
+      delta /= N;
+
+    }
+
+    if (strcmp(method, "graph:") == 0){
+      char* const graph_name = argv[7];
+      cout << "Reading neighborhood graph from file " << graph_name << "..." << endl;
+      pair<double,double> P = compute_delta_from_file(C,graph_name);
+      delta = P.first; lip = P.second;
+    }
+
+    if(VERBOSE)  cout << "lip = " << lip << endl;
+
+    if(VERBOSE)  cout << "delta = " << delta << endl;
+    else cout << delta << endl;
+
+    if(VERBOSE)  cout << "g = " << g << endl;
+    else  cout << g << endl;
+
+    if(!DELTA){
+      if(VERBOSE){
+        if(RESOLUTION)  cout << "r = " << lip*delta/g << endl;
+        else  cout << "r = " << floor(g*(maxf-minf)/(lip*delta*(1-g))) << endl;}
+      else{
+        if(RESOLUTION)  cout << lip*delta/g << endl;
+        else  cout << floor(g*(maxf-minf)/(lip*delta*(1-g))) << endl;}}
+    else{
+      if(VERBOSE){
+        if(RESOLUTION)  cout << "r = " << lip*delta << endl;
+        else  cout << "r = " << floor((maxf-minf)/(lip*delta*(1-g))) << endl;}
+      else{
+        if(RESOLUTION)  cout << lip*delta << endl;
+        else  cout << floor((maxf-minf)/(lip*delta*(1-g))) << endl;}}
+
+    return 0;
+}
diff --git a/src/Nerve_GIC/example/quantile.cpp b/src/Nerve_GIC/example/quantile.cpp
new file mode 100644
index 00000000..14dee5f2
--- /dev/null
+++ b/src/Nerve_GIC/example/quantile.cpp
@@ -0,0 +1,48 @@
+#include "mapper.h"
+
+bool myComp(const pair<double,double> & P1, const pair<double,double> & P2){
+  if(P1.first == P2.first)  return (P1.second <= P2.second);
+  else  return P1.first < P2.first;
+}
+
+int main(int argc, char** argv){
+
+  if(argc <= 3 || argc >= 5){cout << "./quant <file> <F/F^{-1}> <value>" << endl; return 0;}
+
+  char* const file = argv[1];
+  bool token = atoi(argv[2]);
+  double value = atof(argv[3]);
+  vector<double> V;
+  vector<pair<double,double> > Q;
+
+  double v; ifstream input(file); string line;
+  while(getline(input,line)){
+    stringstream stream(line); stream >> v; V.push_back(v);
+  }
+
+  int N = V.size();
+
+  if(!token){
+    int count = 0;
+    for(int i = 0; i < N; i++)
+      if(V[i] <= value)
+        count++;
+    cout << count*1.0/N << endl;
+  }
+  else{
+    for(int i = 0; i < N; i++){
+      int counti = 0;
+      for(int j = 0; j < N; j++)
+        if(V[j] <= V[i])
+          counti++;
+      Q.push_back(pair<double,double>(counti*1.0/N,V[i]));
+    }
+    sort(Q.begin(),Q.end(),myComp);
+    int ind = 0;
+    while(ind < N && Q[ind].first <= value)  ind++;
+    if(ind==N)  cout << Q[ind-1].second << endl;
+    else  cout << Q[ind].second << endl;
+  }
+
+  return 0;
+}
diff --git a/src/Nerve_GIC/example/visu.py b/src/Nerve_GIC/example/visu.py
new file mode 100755
index 00000000..1324e08b
--- /dev/null
+++ b/src/Nerve_GIC/example/visu.py
@@ -0,0 +1,43 @@
+import km
+import numpy as np
+from collections import defaultdict
+
+network = {}
+mapper = km.KeplerMapper(verbose=0)
+data = np.zeros((3,3))
+projected_data = mapper.fit_transform( data, projection="sum", scaler=None )
+
+f = open('mapper.txt','r')
+nodes = defaultdict(list)
+links = defaultdict(list)
+custom = defaultdict(list)
+
+dat = f.readline()
+lens = f.readline()
+param = [float(i) for i in f.readline().split(" ")]
+
+nums = [int(i) for i in f.readline().split(" ")]
+num_nodes = nums[0]
+num_edges = nums[1]
+
+for i in range(0,num_nodes):
+	point = [float(j) for j in f.readline().split(" ")]
+	nodes[  str(int(point[0]))  ] = [  int(point[0]), point[1], int(point[2])  ]
+	links[  str(int(point[0]))  ] = []
+	custom[  int(point[0])  ] = point[1]
+
+m = min([custom[i] for i in range(0,num_nodes)])
+M = max([custom[i] for i in range(0,num_nodes)])
+
+for i in range(0,num_edges):
+	edge = [int(j) for j in f.readline().split(" ")]	
+	links[  str(edge[0])  ].append(  str(edge[1])  )	
+	links[  str(edge[1])  ].append(  str(edge[0])  )
+
+network["nodes"] = nodes
+network["links"] = links
+network["meta"] = lens
+
+mapper.visualize(network, color_function = "Lens average", path_html="mapper_visualization_output.html", title=dat,
+graph_link_distance=30, graph_gravity=0.1, graph_charge=-120, custom_tooltips=custom, width_html=0, 
+height_html=0, show_tooltips=True, show_title=True, show_meta=True, res=param[0],gain=param[1], minimum=m,maximum=M)
\ No newline at end of file
diff --git a/src/Nerve_GIC/include/gudhi/mapper.h b/src/Nerve_GIC/include/gudhi/mapper.h
new file mode 100644
index 00000000..e063b5ae
--- /dev/null
+++ b/src/Nerve_GIC/include/gudhi/mapper.h
@@ -0,0 +1,705 @@
+#include "utils.h"
+
+typedef map<Point, vector<Point> > AdjacencyMatrix; //sparse matrix of adjacencies
+typedef map<int, vector<Point> > ConnectedComp; //a cc is characterized by its ID and the list of its points
+typedef map<int,vector<int> > SparseGraph;
+typedef map<int, pair<double,ConnectedComp> > MapperElements;
+typedef vector<vector<int> > SimplicialComplex;
+
+bool comp(const pair<double, double> & I1, const pair<double, double> & I2){return I1.second < I2.second;}
+
+bool compSC(const vector<int> & V1, const vector<int> & V2){
+  if (V1.size() != V2.size())  return V1.size() < V2.size();
+  else{
+    int tok = 0;
+    while(V1[tok] == V2[tok] && tok < V1.size())  tok++;
+    if(tok != V1.size())  return (V1[tok] < V2[tok]);
+    else return 0;
+  }
+}
+
+class Cover{
+
+  public:
+
+    int res;
+    vector<pair<double, double> > intervals;
+    vector<double> value;
+
+    Cover(){}
+
+    Cover(const double & minf, const double & maxf, const double & resolution, const double & gain, const double & shift, const bool & token){
+      if(!token){
+        double incr = (maxf-minf)/resolution; double x = minf; double alpha = (incr*gain)/(2-2*gain);
+        double y = minf + shift*incr + alpha; pair<double, double> interm(x,y); intervals.push_back(interm);
+        for(int i = 1; i < resolution-1; i++){
+          x = minf + i*incr - alpha;
+          y = minf + (i+1)*incr + alpha;
+          pair<double, double> inter(x,y); intervals.push_back(inter);
+        }
+        x = minf + (resolution-1)*incr - alpha; y = maxf;
+        pair<double, double> interM(x,y); intervals.push_back(interM); res = intervals.size();
+      }
+      else{
+        double x = minf; double y = x + shift*resolution;
+        while(y <= maxf && maxf - (y-gain*resolution) >= resolution){
+          pair<double, double> inter(x,y); intervals.push_back(inter);
+          x = y - gain*resolution;
+          y = x + resolution;
+        }
+        pair<double, double> interM(x,maxf); intervals.push_back(interM); res = intervals.size();
+      }
+    }
+
+    Cover(char* const & name){
+      ifstream input(name); intervals.clear();
+      if(input){
+        string line;
+        while(getline(input,line)){
+          pair<double, double> inter; double x = numeric_limits<double>::max();
+          stringstream stream(line);
+          stream >> x; inter.first = x; stream >> x; inter.second = x;
+          if(x != numeric_limits<double>::max()){assert (inter.first <= inter.second); intervals.push_back(inter);}
+        }
+        res = intervals.size();
+      }
+      else{cout << "  Failed to read file " << name << endl;}
+    }
+
+    void sort_covering(){sort(intervals.begin(),intervals.end(), comp);}
+
+    void proper_value(){
+      value.push_back(0.5*(intervals[0].first+intervals[1].first));
+      for(int i = 1; i < this->res-1; i++)  value.push_back(0.5*(intervals[i-1].second+intervals[i+1].first));
+      value.push_back(0.5*(intervals[res-2].second+intervals[res-1].second));
+    }
+};
+
+bool check_inter(vector<Point> & v1, vector<Point> & v2){
+  vector<Point> v3(v1.size()+v2.size());
+  set_intersection(v1.begin(),v1.end(),v2.begin(),v2.end(),v3.begin());
+  return (v3[0].ID != -1 );
+}
+
+
+
+
+void compute_eccentricity(Cloud & C, char* const & name){
+
+  int num_pts = C.size();
+  double* ecc = new double[num_pts]; for(int i = 0; i < num_pts; i++)  ecc[i] = 0;
+  ifstream input(name, ios::binary);
+  if(input.good()){
+    cout << "  Reading distances..." << endl;
+    double d;
+    for(int i = 0; i < num_pts; i++){
+      for (int j = 0; j < num_pts; j++){
+        //input >> d;
+        input.read((char*) &d,8);
+        if(d >= ecc[i])  ecc[i] = d;
+      }
+    }
+    input.close();
+  }
+  else{cout << "No distance matrix found" << endl; return;}
+
+  for(int i = 0; i < num_pts; i++)  C[i].func = ecc[i];
+  delete [] ecc;
+  return;
+}
+
+pair<double,double> compute_delta_from_file(const Cloud & C, char* const & name){
+
+  ifstream input(name); double maximum = 0; double lip = 0;
+  if(input){
+    string line; int k = 0;
+    while(getline(input,line)){
+      int v = numeric_limits<int>::max(); stringstream stream(line);
+      while(stream >> v){double dis = C[k].EuclideanDistance(C[v]); maximum = max(maximum, dis); lip = max(lip,abs(C[k].func-C[v].func)/dis);}
+      k++;
+    }
+    return pair<double,double>(maximum,lip);
+  }
+  else{cout << "Failed to read file " << name << endl; return pair<double,double>(0,0);}
+
+}
+
+AdjacencyMatrix build_neighborhood_graph_from_file(const Cloud & C, char* const & name){
+
+  int nb = C.size(); AdjacencyMatrix G; G.clear(); vector<Point> adj; adj.clear();
+  for(int i = 0; i < nb; i++)
+    G.insert(pair<Point, vector<Point> >(C[i],adj));
+  ifstream input(name);
+  if(input){
+    string line; int num_edges = 0; int k = 0;
+    while(getline(input,line)){
+      vector<Point> ad; ad.clear(); int v = numeric_limits<int>::max(); stringstream stream(line);
+      while(stream >> v)  ad.push_back(C[v]);
+      if(v != numeric_limits<int>::max()){G[C[k]] = ad; k++; num_edges += ad.size();}
+    }
+    cout << "  " << 100*((double) num_edges/2)/((double) nb*(nb-1)/2) << "% of pairs selected." << endl;
+    return G;
+  }
+  else{cout << "Failed to read file " << name << endl; return G;}
+
+}
+
+AdjacencyMatrix build_neighborhood_graph_from_scratch_with_percentage(const Cloud & C, const double & delta, char* const & name, const bool & VNE){
+
+  ifstream input(name, ios::out | ios::binary);
+
+  if(input.good()){
+    cout << "  Reading distances..." << endl;
+    int nb = C.size(); AdjacencyMatrix G; G.clear(); vector<Point> adj;
+    double d; vector<vector<double> > dist; dist.clear();
+
+    double m = 0;
+    for(int i = 0; i < nb; i++){
+      vector<double> dis; dis.clear();
+      for (int j = 0; j < nb; j++){
+        input.read((char*) &d,8); dis.push_back(d);
+        if(m <= d)  m = d;
+      }
+      dist.push_back(dis);
+    }
+    input.close();
+
+    cout << "  Done." << endl << "  Computing neighborhood graph...  ";
+
+    cout.flush(); int num_edges = 0;
+    for(int i = 0; i < nb; i++){
+      adj.clear();
+      for (int j = 0; j < nb; j++)
+        if(dist[i][j] <= delta*m){adj.push_back(C[j]); num_edges++;}
+      G.insert(pair<Point, vector<Point> >(C[i],adj));
+    }
+    cout << 100*((double) num_edges/2)/((double) nb*(nb-1)/2) << "% of pairs selected." << endl;
+    return G;
+  }
+  else{
+      cout << "  Computing distances..." << endl; vector<double> V; double m = 0;
+      input.close(); ofstream output(name, ios::out | ios::binary);
+      int nb = C.size(); AdjacencyMatrix G; G.clear(); vector<Point> adj;
+      double d; vector<vector<double> > dist; dist.clear();
+
+      if(VNE){
+        int dim = C[0].coord.size(); int N = C.size();
+        for(int i = 0; i < dim; i++){
+          double meani = 0;
+          for (int j = 0; j < nb; j++)  meani += C[j].coord[i]/N;
+          double vari = 0;
+          for (int j = 0; j < nb; j++)  vari += pow(C[j].coord[i]-meani,2)/N;
+          V.push_back(vari);
+        }
+      }
+
+      for(int i = 0; i < nb; i++){
+        if( (int) floor( 100*((double) i)/((double) nb)+1 ) %10 == 0  ){cout << "\r" << floor( 100*((double) i)/((double) nb) +1) << "%" << flush;}
+        vector<double> dis; dis.clear();
+        for (int j = 0; j < nb; j++){
+          if(!VNE){d = C[i].EuclideanDistance(C[j]); dis.push_back(d); if(m <= d)  m = d;}
+          else{d = C[i].VarianceNormalizedEuclideanDistance(C[j],V); dis.push_back(d); if(m <= d)  m = d;}
+          output.write((char*) &d,8);
+        }
+        dist.push_back(dis);
+      }
+      output.close();
+
+    cout << endl << "  Done." << endl << "  Computing neighborhood graph...  ";
+
+    cout.flush(); int num_edges = 0;
+    for(int i = 0; i < nb; i++){
+      adj.clear();
+      for (int j = 0; j < nb; j++)
+        if(dist[i][j] <= delta*m){adj.push_back(C[j]); num_edges++;}
+      G.insert(pair<Point, vector<Point> >(C[i],adj));
+    }
+    cout << 100*((double) num_edges)/pow(nb,2) << "% of pairs selected." << endl;
+    return G;
+  }
+}
+
+
+
+AdjacencyMatrix build_neighborhood_graph_from_scratch_with_parameter(const Cloud & C, const double & delta, char* const & name, const bool & VNE){
+
+    ifstream input(name, ios::out | ios::binary);
+
+    if(input.good()){
+      cout << "  Reading distances..." << endl;
+      int nb = C.size(); AdjacencyMatrix G; G.clear(); vector<Point> adj;
+      double d; vector<vector<double> > dist; dist.clear();
+
+      for(int i = 0; i < nb; i++){
+        vector<double> dis; dis.clear();
+        for (int j = 0; j < nb; j++){
+          input.read((char*) &d,8); dis.push_back(d);
+        }
+        dist.push_back(dis);
+      }
+      input.close();
+
+      cout << "  Done." << endl << "  Computing neighborhood graph...  ";
+
+      cout.flush(); int num_edges = 0;
+      for(int i = 0; i < nb; i++){
+        adj.clear();
+        for (int j = 0; j < nb; j++)
+          if(dist[i][j] <= delta && j != i){adj.push_back(C[j]); num_edges++;}
+        G.insert(pair<Point, vector<Point> >(C[i],adj));
+      }
+      cout << 100*((double) num_edges/2)/((double) nb*(nb-1)/2) << "% of pairs selected." << endl;
+      return G;
+    }
+    else{
+      cout << "  Computing distances..." << endl; vector<double> V;
+      input.close(); ofstream output(name, ios::out | ios::binary);
+      int nb = C.size(); AdjacencyMatrix G; G.clear(); vector<Point> adj;
+      double d; vector<vector<double> > dist; dist.clear();
+
+      if(VNE){
+        int dim = C[0].coord.size(); int N = C.size();
+        for(int i = 0; i < dim; i++){
+          double meani = 0;
+          for (int j = 0; j < nb; j++)  meani += C[j].coord[i]/N;
+          double vari = 0;
+          for (int j = 0; j < nb; j++)  vari += pow(C[j].coord[i]-meani,2)/N;
+          V.push_back(vari);
+        }
+      }
+
+      for(int i = 0; i < nb; i++){
+        if( (int) floor( 100*((double) i)/((double) nb)+1 ) %10 == 0  ){cout << "\r" << floor( 100*((double) i)/((double) nb) +1) << "%" << flush;}
+        vector<double> dis; dis.clear();
+        for (int j = 0; j < nb; j++){
+          if(!VNE){d = C[i].EuclideanDistance(C[j]); dis.push_back(d);}
+          else{d = C[i].VarianceNormalizedEuclideanDistance(C[j],V); dis.push_back(d);}
+          output.write((char*) &d,8);
+        }
+        dist.push_back(dis);
+      }
+      output.close();
+
+      cout << endl << "  Done." << endl << "  Computing neighborhood graph...  ";
+
+      cout.flush(); int num_edges = 0;
+      for(int i = 0; i < nb; i++){
+        adj.clear();
+        for (int j = 0; j < nb; j++)
+          if(dist[i][j] <= delta && j != i){adj.push_back(C[j]); num_edges++;}
+        G.insert(pair<Point, vector<Point> >(C[i],adj));
+      }
+      cout << 100*((double) num_edges)/pow(nb,2) << "% of pairs selected." << endl;
+      return G;
+    }
+}
+
+
+
+
+void dfs(AdjacencyMatrix & G, const Point & p, vector<Point> & cc, map<Point,bool> & visit){
+  cc.push_back(p);
+  visit[p] = true; int neighb = G[p].size();
+  for (int i = 0; i < neighb; i++)
+    if (  visit.find(G[p][i]) != visit.end() )
+      if(  !(visit[G[p][i]])  )
+        dfs(G,G[p][i],cc,visit);
+}
+
+ConnectedComp count_cc(AdjacencyMatrix & G, int & id, map<Point,vector<int> > & clusters){
+  map<Point,bool> visit;
+  for(AdjacencyMatrix::iterator it = G.begin(); it != G.end(); it++)
+    visit.insert(pair<Point,bool>(it->first, false));
+  ConnectedComp CC; CC.clear();
+  if (!(G.empty())){
+    for(AdjacencyMatrix::iterator it = G.begin(); it != G.end(); it++){
+      if (  !(visit[it->first])  ){
+        vector<Point> cc; cc.clear();
+        dfs(G,it->first,cc,visit); int cci = cc.size();
+        for(int i = 0; i < cci; i++)  clusters[cc[i]].push_back(id);
+        CC.insert(pair<int, vector<Point> >(id++,cc));
+      }
+    }
+  }
+  return CC;
+}
+
+MapperElements MapperElts(AdjacencyMatrix & G, const Cover & I, map<Point,vector<int> > & clusters, vector<double> & col){
+
+  map<int, pair<double,ConnectedComp> > mapper_elts;
+  AdjacencyMatrix::iterator pos = G.begin();
+  int id = 0;
+
+  for(int i = 0; i < I.res; i++){
+
+    AdjacencyMatrix prop; prop.clear();
+    pair<double, double> inter1 = I.intervals[i];
+    AdjacencyMatrix::iterator tmp = pos;
+
+    if(i != I.res-1){
+      if(i != 0){
+        pair<double, double> inter3 = I.intervals[i-1];
+        while(tmp->first.func < inter3.second && tmp != G.end()){
+          col[tmp->first.ID] = 0;
+          prop.insert(*tmp);
+          tmp++;
+        }
+      }
+      pair<double, double> inter2 = I.intervals[i+1];
+      while(tmp->first.func < inter2.first && tmp != G.end()){
+        col[tmp->first.ID] = i+1;
+        prop.insert(*tmp);
+        tmp++;
+      }
+      pos = tmp;
+      while(tmp->first.func < inter1.second && tmp != G.end()){
+        prop.insert(*tmp);
+        tmp++;
+      }
+
+    }
+    else{
+      pair<double, double> inter3 = I.intervals[i-1];
+      while(tmp->first.func < inter3.second && tmp != G.end()){
+        col[tmp->first.ID] = 0;
+        prop.insert(*tmp);
+        tmp++;
+      }
+      while(tmp != G.end()){
+        col[tmp->first.ID] = i+1;
+        prop.insert(*tmp);
+        tmp++;
+      }
+
+    }
+    ConnectedComp CC = count_cc(prop,id,clusters);
+    mapper_elts.insert(pair<int, pair<double,ConnectedComp> >(i, pair<double,ConnectedComp>(I.value[i],CC)));
+
+  }
+
+  return mapper_elts;
+}
+
+void find_all_simplices(vector<vector<int> > & lc, SimplicialComplex & SC, int & tok, vector<int> & simpl){
+  int num_nodes = lc.size();
+  if(tok == num_nodes-1){
+    int num_clus = lc[tok].size();
+    for(int i = 0; i < num_clus; i++){
+      vector<int> simplex = simpl; simplex.push_back(lc[tok][i]);
+      sort(simplex.begin(), simplex.end()); vector<int>::iterator iter = unique(simplex.begin(),simplex.end()); simplex.resize(distance(simplex.begin(),iter));
+      SC.push_back(simplex); int dimension = simplex.size();
+      if(dimension >= 2){for(int j = 0; j < dimension; j++){vector<int> face = simplex; face.erase(face.begin()+j); SC.push_back(face);}}
+    }
+  }
+  else{
+    int num_clus = lc[tok].size();
+    for(int i = 0; i < num_clus; i++){
+      vector<int> simplex = simpl; simplex.push_back(lc[tok][i]);
+      find_all_simplices(lc, SC, ++tok, simplex);
+    }
+  }
+}
+
+void bron_kerbosch_II(vector<Point> R, vector<Point> P, vector<Point> X, AdjacencyMatrix & G, vector<vector<Point> > & cliques){
+
+  if(P.size() == 0 && X.size() == 0){cliques.push_back(R); /*cout << "clique "; for(int i = 0; i < R.size(); i++)  cout << R[i].ID+1 << "  "; cout << endl;*/}
+  else{
+    if(P.size() != 0){
+      Point pivot = P[0]; vector<Point> neighbors_pivot = G[pivot];
+
+      sort(P.begin(),P.end()); //cout << "P = "; for(int i = 0; i < P.size(); i++)  cout << P[i].ID+1 << "  "; cout << endl;
+      sort(neighbors_pivot.begin(),neighbors_pivot.end()); //cout << "pivot = " << pivot.ID+1 << endl;
+
+      vector<Point> list(neighbors_pivot.size()+P.size()); vector<Point>::iterator it;
+      it = set_difference(P.begin(),P.end(),neighbors_pivot.begin(),neighbors_pivot.end(),list.begin());
+      list.resize(distance(list.begin(),it)); int nb = list.size();
+      //cout << "N(pivot) = "; for(int i = 0; i < neighbors_pivot.size(); i++)  cout << G[pivot][i].ID+1 << "  "; cout << endl;
+      //cout << "P minus N(pivot) = "; for(int i = 0; i < list.size(); i++)  cout << list[i].ID+1 << "  "; cout << endl;
+
+      for(int i = 0; i < nb; i++){
+
+        Point v = list[i]; //cout << "v = " << v.ID+1 << endl;
+        vector<Point> neighbors_v = G[v]; sort(neighbors_v.begin(),neighbors_v.end());
+        vector<Point> newR = R; newR.push_back(v);
+
+        //cout << "P = "; for(int j = 0; j < P.size(); j++)  cout << P[j].ID+1 << "  "; cout << endl;
+        //cout << "N(v) = "; for(int j = 0; j < neighbors_v.size(); j++)  cout << neighbors_v[j].ID+1 << "  "; cout << endl;
+
+        vector<Point> newP(P.size()+neighbors_v.size()); it = set_intersection(P.begin(),P.end(),neighbors_v.begin(),neighbors_v.end(),newP.begin());
+        newP.resize(it-newP.begin()); //cout << "new P = "; for(int j = 0; j < newP.size(); j++)  cout << newP[j].ID+1 << "  "; cout << endl;
+
+        vector<Point> newX(X.size()+neighbors_v.size()); it = set_intersection(X.begin(),X.end(),neighbors_v.begin(),neighbors_v.end(),newX.begin());
+        newX.resize(it-newX.begin()); //cout << "new X = "; for(int j = 0; j < newX.size(); j++)  cout << newX[j].ID+1 << "  "; cout << endl;
+
+        bron_kerbosch_II(newR,newP,newX,G,cliques);
+        it = find(P.begin(),P.end(),v); P.erase(it);
+        X.push_back(v);
+
+      }
+    }
+  }
+}
+
+SimplicialComplex compute_GIC(AdjacencyMatrix & G, map<Point,vector<int> > & clusters){
+  SimplicialComplex SC; AdjacencyMatrix H;
+  vector<Point> R, P, X; R.clear(); X.clear();
+  for(AdjacencyMatrix::iterator it = G.begin(); it!=G.end(); it++){
+    int numneigh = it->second.size(); int numclus = clusters[it->first].size();
+    if(numclus > 1){H.insert(*it); P.push_back(it->first);}
+    else{
+      int currentcluster = clusters[it->first][0]; vector<Point> vectneigh;
+      for(int i = 0; i < numneigh; i++){
+        if(clusters[it->second[i]].size() > 1 || clusters[it->second[i]][0] != currentcluster)
+          vectneigh.push_back(it->second[i]);
+      }
+      if(vectneigh.size() > 0){
+        H.insert(pair<Point,vector<Point> >(it->first, vectneigh));
+        P.push_back(it->first);
+      }
+    }
+    vector<int> vertex(1); for(int i = 0; i < numclus; i++){vertex[0] = clusters[it->first][i]; SC.push_back(vertex);}
+  }
+  vector<vector<Point> > cliques; cout << "  Computing maximal cliques..." << endl;
+  bron_kerbosch_II(R,P,X,H,cliques); cout << "  Done" << endl;
+  int num_cliques = cliques.size();
+  for(int i = 0; i < num_cliques; i++){
+    int num_nodes = cliques[i].size(); vector<vector<int> > list_clusters(num_nodes);
+    for(int j = 0; j < num_nodes; j++)  list_clusters[j] = clusters[cliques[i][j]];
+    vector<int> simpl; int start = 0; find_all_simplices(list_clusters,SC,start,simpl);
+  }
+  sort(SC.begin(), SC.end(), compSC); SimplicialComplex::iterator iter = unique(SC.begin(), SC.end()); SC.resize(iter-SC.begin());
+  return SC;
+}
+
+void add_simplices(vector<int> & clus, SimplicialComplex & SC){
+  int dimension = clus.size();
+  if(dimension == 1)  SC.push_back(clus);
+  else{
+    SC.push_back(clus);
+    for(int i = 0; i < dimension; i++){
+      vector<int> cl = clus; cl.erase(cl.begin()+i);
+      add_simplices(cl,SC);
+    }
+  }
+}
+
+SimplicialComplex compute_Nerve(vector<Point> & C, map<Point,vector<int> > & clusters){
+  SimplicialComplex SC; int num_pts = C.size();
+  for(int i = 0; i < num_pts; i++){
+    vector<int> clus = clusters[C[i]];
+    add_simplices(clus,SC);
+  }
+  sort(SC.begin(), SC.end(), compSC); SimplicialComplex::iterator iter = unique(SC.begin(), SC.end()); SC.resize(iter-SC.begin());
+  return SC;
+}
+
+SparseGraph build_mapper_graph(MapperElements & M, map<int, double> & colors, map<int,int> & num){
+
+  SparseGraph MG; MG.clear(); vector<int> adj;
+  int nb_elts = M.size(); adj.clear();
+
+  #pragma omp parallel for
+  for (int i = 0; i < nb_elts; i++){
+    double c = M[i].first;
+    for (ConnectedComp::iterator it = M[i].second.begin(); it != M[i].second.end(); it++){
+      int p = it->first;
+      #pragma omp critical
+      {
+        MG.insert(pair<int, vector<int> >(p,adj));
+        int inside_nb = it->second.size();
+        colors.insert(pair<int,double>(p,c)); num.insert(pair<int,int>(p,inside_nb));
+      }
+    }
+  }
+
+
+  for(int i = 0; i < nb_elts-1; i++){
+
+    if(i==0){
+      for (ConnectedComp::iterator it = M[i].second.begin(); it != M[i].second.end(); it++)
+        sort(it->second.begin(),it->second.end());
+    }
+
+    for (ConnectedComp::iterator iit = M[i+1].second.begin(); iit != M[i+1].second.end(); iit++)
+      sort(iit->second.begin(),iit->second.end());
+
+    for (ConnectedComp::iterator it = M[i].second.begin(); it != M[i].second.end(); it++)
+      for (ConnectedComp::iterator iit = M[i+1].second.begin(); iit != M[i+1].second.end(); iit++)
+        if(check_inter(it->second, iit->second))
+          MG[it->first].push_back(iit->first);
+  }
+
+  return MG;
+}
+
+SparseGraph build_mapperDelta_graph(MapperElements & M, AdjacencyMatrix & G, map<Point,vector<int> > & clusters, \
+                                    map<int, double> & colors, map<int,int> & num){
+
+  SparseGraph MG; MG.clear(); vector<int> adj;
+  int nb_elts = M.size(); adj.clear();
+
+  #pragma omp parallel for
+  for (int i = 0; i < nb_elts; i++){
+    double c = M[i].first;
+    for (ConnectedComp::iterator it = M[i].second.begin(); it != M[i].second.end(); it++){
+      int p = it->first;
+      #pragma omp critical
+      {
+        MG.insert(pair<int, vector<int> >(p,adj));
+        int inside_nb = it->second.size();
+        colors.insert(pair<int,double>(p,c)); num.insert(pair<int,int>(p,inside_nb));
+      }
+    }
+  }
+
+  #pragma omp parallel for
+  for(int i = 0; i < nb_elts-1; i++){
+    double vali = M[i+1].first;
+    for (ConnectedComp::iterator it = M[i].second.begin(); it != M[i].second.end(); it++){
+      int sit = it->second.size(); int idit = it->first;
+      for(int j = 0; j < sit; j++){
+        Point p = it->second[j]; int neighp = G[p].size(); vector<int> clusp = clusters[p]; int numclus = clusp.size();
+        for(int k = 0; k < numclus; k++)  if(clusp[k] > idit)  MG[idit].push_back(clusp[k]);
+        for(int k = 0; k < neighp; k++){
+          Point q = G[p][k]; int clusq = clusters[q].size();
+          for(int l = 0; l < clusq; l++)
+            if(clusters[q][l] > idit && colors[clusters[q][l]] == vali){
+            #pragma omp critical
+            {MG[idit].push_back(clusters[q][l]);}
+            }
+        }
+      }
+      sort(MG[idit].begin(),MG[idit].end()); vector<int>::iterator iter = unique(MG[idit].begin(),MG[idit].end());
+      MG[idit].resize(distance(MG[idit].begin(),iter));
+    }
+  }
+
+  return MG;
+}
+
+/*
+vector<pair<int,int> > check_intersection_crossing(const vector<ConnectedComp> & M, AdjacencyMatrix & G){
+  int nb_elts = M.size(); vector<pair<int,int> > pairs; pairs.clear();
+  for (int i = 0; i < nb_elts-1; i+=2){
+    ConnectedComp CC1 = M[i]; ConnectedComp CC2 = M[i+2];
+    for(ConnectedComp::iterator it = CC1.begin(); it != CC1.end(); it++){
+      int cc1_nb = it->second.size();
+      for(ConnectedComp::iterator iit = CC2.begin(); iit != CC2.end(); iit++){
+        int cc2_nb = iit->second.size(); int i = 0; bool point_in_inter = false; int count_interval_cross = 0;
+        //cout << "  Checking nodes " << it->first << " and " << iit->first << "..." << endl;
+        while(i < cc1_nb && !point_in_inter){
+          Point p1 = it->second[i]; int j = 0;
+          while(j < cc2_nb && !point_in_inter){
+            Point p2 = iit->second[j];
+            if(p1 == p2){point_in_inter = true;}
+            if(  find(G[p1].begin(),G[p1].end(),p2)!=G[p1].end()  ){count_interval_cross++;}
+            j++;
+          }
+          i++;
+        }
+        if(!point_in_inter && count_interval_cross > 0){
+          cout << "  Warning: nodes " << it->first << " and " << iit->first << " separated whereas they should not: " << count_interval_cross << " intersection crossings" << endl;
+          pairs.push_back(pair<int,int>(it->first,iit->first));
+        }
+      }
+    }
+  }
+  return pairs;
+}
+*/
+
+void plotNG(const int & dim, AdjacencyMatrix & G){
+  char neigh[100]; int ned = 0; for (AdjacencyMatrix::iterator it = G.begin(); it != G.end(); it++){ned += it->second.size();}
+  sprintf(neigh, "NG.vect");
+  ofstream output(neigh);
+  output << "VECT" << endl;
+  output << ned + G.size() << " " << 2*ned + G.size() << " " << ned + G.size() << endl;
+  for(int i = 0; i < ned; i++)  output << "2 ";
+  for(int i = 0; i < G.size(); i++)  output << "1 ";
+  output << endl;
+  for(int i = 0; i < ned + G.size(); i++)  output << "1 ";
+  output << endl;
+  for (AdjacencyMatrix::iterator it = G.begin(); it != G.end(); it++){
+    for(int j = 0; j < it->second.size(); j++){
+      for(int k = 0; k < dim; k++)  output << it->first.coord[k] << " ";
+      if(dim == 1)  output << "0 0  ";
+      if(dim == 2)  output << "0  ";
+      if(dim == 3)  output << "  ";
+      for(int k = 0; k < dim; k++)  output << it->second[j].coord[k] << " ";
+      if(dim == 1)  output << "0 0" << endl;
+      if(dim == 2)  output << "0" << endl;
+      if(dim == 3)  output << endl;
+    }
+  }
+  for (AdjacencyMatrix::iterator it = G.begin(); it != G.end(); it++){
+    for(int k = 0; k < dim; k++)  output << it->first.coord[k] << " ";
+    if(dim == 1)  output << "0 0" << endl;
+    if(dim == 2)  output << "0" << endl;
+    if(dim == 3)  output << endl;
+  }
+  //for(int i = 0; i < ned + G.size(); i++){
+  int plot_coord = 0;
+  double c; double maxc = G.begin()->first.coord[2]; double minc = G.begin()->first.coord[plot_coord];
+  if(dim == 3){
+    for (AdjacencyMatrix::iterator it = G.begin(); it != G.end(); it++){
+      if(it->first.coord[plot_coord] >= maxc)  maxc = it->first.coord[plot_coord];
+      if(it->first.coord[plot_coord] <= minc)  minc = it->first.coord[plot_coord];
+    }
+  }
+  for (AdjacencyMatrix::iterator it = G.begin(); it != G.end(); it++){
+    if(dim == 3)  c = (it->first.coord[plot_coord]-minc)/(maxc-minc);
+    else  c = 0;
+    for(int j = 0; j < it->second.size(); j++){
+      if(c<=0.5)  output << 1-2*c << " " << 2*c << " " << 0 << " 1" << endl;
+      else  output << 0 << " " << -2*c+2 << " " << 2*c-1 << " 1" << endl;
+    }
+  }
+  for (AdjacencyMatrix::iterator it = G.begin(); it != G.end(); it++){
+    if(dim == 3)  c = (it->first.coord[plot_coord]-minc)/(maxc-minc);
+    else  c = 0;
+    if(c<=0.5)  output << 1-2*c << " " << 2*c << " " << 0 << " 1" << endl;
+    else  output << 0 << " " << -2*c+2 << " " << 2*c-1 << " 1" << endl;
+  }
+  //   output << "0 0 1 1" << endl;
+}
+
+void plotCover(const int & dim, AdjacencyMatrix & G, const vector<double> & col, const Cover & I){
+  int number = G.size();
+  char neigh[100]; sprintf(neigh, "Cover.vect"); ofstream output(neigh);
+  output << "VECT" << endl;
+  output << number << " " << number << " " << number << endl;
+  for(int i = 0; i < number; i++)  output << "1 ";
+  output << endl;
+  for(int i = 0; i < number; i++)  output << "1 ";
+  output << endl;
+  for (AdjacencyMatrix::iterator it = G.begin(); it != G.end(); it++){
+    for(int k = 0; k < dim; k++)  output << it->first.coord[k] << " ";
+    if(dim == 1)  output << "0 0" << endl;
+    if(dim == 2)  output << "0" << endl;
+    if(dim == 3)  output << endl;
+  }
+  int co = 1;
+  for (AdjacencyMatrix::iterator it = G.begin(); it != G.end(); it++){
+    double c = col[it->first.ID]/I.res;
+    if(co != G.size()){
+      if(c<=0.5){
+        if(c==0)  output << 0 << " " << 0 << " " << 0 << " 1" << endl;
+        else{
+          if(c==1.0/I.res)  output << 1 << " " << 0 << " " << 0 << " 1" << endl;
+          else  output << 1-2*c << " " << 2*c << " " << 0 << " 1" << endl;
+        }
+      }
+      else  output << 0 << " " << -2*c+2 << " " << 2*c-1 << " 1" << endl;
+      co++;}
+    else{
+      if(c<=0.5){
+        if(c==0)  output << 0 << " " << 0 << " " << 0 << " 1" << endl;
+        else{
+          if(c==1.0/I.res)  output << 1 << " " << 0 << " " << 0 << " 1" << endl;
+          else  output << 1-2*c << " " << 2*c << " " << 0 << " 1" << endl;
+        }
+      }
+      else  output << 0 << " " << -2*c+2 << " " << 2*c-1 << " 1" << endl;
+    }
+  }
+}
diff --git a/src/Nerve_GIC/include/gudhi/utils.h b/src/Nerve_GIC/include/gudhi/utils.h
new file mode 100644
index 00000000..55b54991
--- /dev/null
+++ b/src/Nerve_GIC/include/gudhi/utils.h
@@ -0,0 +1,86 @@
+#include <cstdlib>
+#include <string.h>
+#include <stdio.h>
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <vector>
+#include <algorithm>
+#include <set>
+#include <map>
+#include <limits>
+#include <assert.h>
+#include <math.h>
+#include <iterator>
+#include <omp.h>
+#include <list>
+#include <random>
+
+using namespace std;
+
+class Point{
+
+  public:
+
+    int ID;
+    double func;
+    vector<double> coord;
+
+    Point(){ ID = -1;}
+    Point(const int & _ID){ID = _ID;} // constructor 1
+    Point(const int & _ID, const vector<double> & _coord){ID = _ID; coord = _coord;} // constructor 2
+    bool operator<(const Point & p) const {if(func != p.func){return func < p.func;}else{return ID < p.ID;}} // comparator
+    bool operator==(const Point & p) const {return ID == p.ID;}
+    double EuclideanDistance(const Point & p) const {
+      //cout << this->coord.size() << " " << p.coord.size() << endl;
+      assert (coord.size() == p.coord.size());
+      double x = 0; int dim = p.coord.size(); for (int i = 0; i < dim; i++){x+=(coord[i]-p.coord[i])*(coord[i]-p.coord[i]);}
+      return sqrt(x);
+    }
+    double VarianceNormalizedEuclideanDistance(const Point & p, const vector<double> & V) const {
+      //cout << this->coord.size() << " " << p.coord.size() << endl;
+      assert (coord.size() == p.coord.size());
+      double x = 0; int dim = p.coord.size(); for (int i = 0; i < dim; i++){x+=(coord[i]-p.coord[i])*(coord[i]-p.coord[i])/V[i];}
+      return sqrt(x);
+    }
+
+};
+
+typedef vector<Point> Cloud;
+
+Cloud read_cloud(char* const & name){
+  int dim; Cloud C; C.clear();
+  ifstream input(name);
+  if(input){
+    string line; vector<double> coord; int ID = 0;
+    while(getline(input,line)){
+      coord.clear(); double x = numeric_limits<double>::max();
+      stringstream stream(line);
+      while(stream >> x){coord.push_back(x);}
+      Point p = Point(ID, coord);
+      if(x != numeric_limits<double>::max()){C.push_back(p); ID++;}
+    }
+  }
+  else{cout << "Failed to read file " << name << endl; return C;}
+  return C;
+}
+
+void read_function_from_file(char* const & name, Cloud & C){
+  int num_pts = C.size(); double x;
+  ifstream input(name);
+  if(input){
+    for(int i = 0; i < num_pts; i++)
+      input >> C[i].func;
+  }
+  else{cout << "Failed to read file " << name << endl; return;}
+  return;
+}
+
+void read_coordinate(const int & number, Cloud & C){
+  int num_pts = C.size();
+  for(int i = 0; i < num_pts; i++)
+    C[i].func = C[i].coord[number];
+  return;
+}
+
+
-- 
cgit v1.2.3