From 69e3499fb9ecae20a1b9904cac83bbc605972c50 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Thu, 16 Nov 2017 20:52:37 +0000 Subject: Apply google code conventions for cpplint git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/Nerve_GIC@2898 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 2f82ddf869dc95ec33c6cea825ee24cd1a849c4b --- src/Nerve_GIC/include/gudhi/GIC.h | 1939 ++++++++++++++++++++----------------- 1 file changed, 1066 insertions(+), 873 deletions(-) (limited to 'src/Nerve_GIC/include') diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h index ca8727a0..f2a01e97 100644 --- a/src/Nerve_GIC/include/gudhi/GIC.h +++ b/src/Nerve_GIC/include/gudhi/GIC.h @@ -35,8 +35,8 @@ #include #include #include -#include // for numeric_limits -#include // for pair<> +#include // for numeric_limits +#include // for pair<> #include // for std::max #include #include @@ -49,7 +49,6 @@ using Simplex_tree = Gudhi::Simplex_tree<>; using Filtration_value = Simplex_tree::Filtration_value; using Rips_complex = Gudhi::rips_complex::Rips_complex; - /** * \class Cover_complex * \brief Cover complex data structure. @@ -70,906 +69,1100 @@ using Rips_complex = Gudhi::rips_complex::Rips_complex; * correspond to the simplices of the GIC. * */ -template +template class Cover_complex { - private: - //Graph_induced_complex(std::map fun){func = fun;} - bool verbose = false; // whether to display information. - std::vector point_cloud; - std::vector > one_skeleton; - typedef int Cover_t; // elements of cover C are indexed by integers. - std::vector > simplices; - std::map > cover; - std::map > cover_back; - int maximal_dim; // maximal dimension of output simplicial complex. - int data_dimension; // dimension of input data. - int n; // number of points. - std::map cover_fct; // integer-valued function that allows to state if two elements of the cover are consecutive or not. - std::map > cover_color; // size and coloring of the vertices of the output simplicial complex. - Simplex_tree st; - - std::map > adjacency_matrix; - std::vector > distances; - - int resolution_int = -1; - double resolution_double = -1; - double gain = -1; - double rate_constant = 10; // Constant in the subsampling. - double rate_power = 0.001; // Power in the subsampling. - int mask = 0; // Ignore nodes containing less than mask points. - - std::map func; - std::map func_color; - std::vector voronoi_subsamples; - std::string cover_name; - std::string point_cloud_name; - std::string color_name; - std::string type; // Nerve or GIC - bool functional_cover = false; // whether we use a cover with preimages of a function or not - - // Point comparator - struct Less{ - Less(std::map func){Fct = func;} - std::map Fct; - bool operator()(int a, int b){if(Fct[a] == Fct[b]) return a < b; else return Fct[a] < Fct[b];} - }; - - // DFS + // Graph_induced_complex(std::map fun){func = fun;} + bool verbose = false; // whether to display information. + std::vector point_cloud; + std::vector > one_skeleton; + typedef int Cover_t; // elements of cover C are indexed by integers. + std::vector > simplices; + std::map > cover; + std::map > cover_back; + int maximal_dim; // maximal dimension of output simplicial complex. + int data_dimension; // dimension of input data. + int n; // number of points. + std::map + cover_fct; // integer-valued function that allows to state if two elements of the cover are consecutive or not. + std::map > + cover_color; // size and coloring of the vertices of the output simplicial complex. + Simplex_tree st; + + std::map > adjacency_matrix; + std::vector > distances; + + int resolution_int = -1; + double resolution_double = -1; + double gain = -1; + double rate_constant = 10; // Constant in the subsampling. + double rate_power = 0.001; // Power in the subsampling. + int mask = 0; // Ignore nodes containing less than mask points. + + std::map func; + std::map func_color; + std::vector voronoi_subsamples; + std::string cover_name; + std::string point_cloud_name; + std::string color_name; + std::string type; // Nerve or GIC + bool functional_cover = false; // whether we use a cover with preimages of a function or not + + // Point comparator + struct Less { + Less(std::map func) { Fct = func; } + std::map Fct; + bool operator()(int a, int b) { + if (Fct[a] == Fct[b]) + return a < b; + else + return Fct[a] < Fct[b]; + } + }; + + // DFS private: - void dfs(std::map >& G, int p, std::vector& cc, std::map& 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); - } - - // Find random number in [0,1]. - double GetUniform(){ - static std::default_random_engine re; - static std::uniform_real_distribution Dist(0,1); - return Dist(re); - } - - // Subsample points. - void SampleWithoutReplacement(int populationSize, int sampleSize, std::vector & samples){ - int t = 0; int m = 0; double u; - while (m < sampleSize){ - u = GetUniform(); - if ( (populationSize - t)*u >= sampleSize - m ) - t++; - else{samples[m] = t; t++; m++;} - } - } + void dfs(std::map >& G, int p, std::vector& cc, std::map& 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); + } + + // Find random number in [0,1]. + double GetUniform() { + static std::default_random_engine re; + static std::uniform_real_distribution Dist(0, 1); + return Dist(re); + } + + // Subsample points. + void SampleWithoutReplacement(int populationSize, int sampleSize, std::vector& samples) { + int t = 0; + int m = 0; + double u; + while (m < sampleSize) { + u = GetUniform(); + if ((populationSize - t) * u >= sampleSize - m) + t++; + else { + samples[m] = t; + t++; + m++; + } + } + } private: - void fill_adjacency_matrix_from_st(){ - std::vector empty; - for(int i = 0; i < n; i++) adjacency_matrix[i] = empty; - for (auto simplex : st.complex_simplex_range()) { - if(st.dimension(simplex) == 1){ - std::vector vertices; - for(auto vertex : st.simplex_vertex_range(simplex)) vertices.push_back(vertex); - adjacency_matrix[vertices[0]].push_back(vertices[1]); adjacency_matrix[vertices[1]].push_back(vertices[0]); - } - } - } - - public: - /** \brief Specifies whether the type of the output simplicial complex. - * - * @param[in] t string (either "GIC" or "Nerve"). - * - */ - void set_type(const std::string & t){type = t;} + void fill_adjacency_matrix_from_st() { + std::vector empty; + for (int i = 0; i < n; i++) adjacency_matrix[i] = empty; + for (auto simplex : st.complex_simplex_range()) { + if (st.dimension(simplex) == 1) { + std::vector vertices; + for (auto vertex : st.simplex_vertex_range(simplex)) vertices.push_back(vertex); + adjacency_matrix[vertices[0]].push_back(vertices[1]); + adjacency_matrix[vertices[1]].push_back(vertices[0]); + } + } + } public: - /** \brief Specifies whether the program should display information or not. - * - * @param[in] verb boolean (true = display info, false = do not display info). - * - */ - void set_verbose(bool verb = false){verbose = verb;} + /** \brief Specifies whether the type of the output simplicial complex. + * + * @param[in] t string (either "GIC" or "Nerve"). + * + */ + void set_type(const std::string& t) { type = t; } public: - /** \brief Sets the constants used to subsample the data set. These constants are - * explained in \cite Carriere17c. - * - * @param[in] constant double. - * @param[in] power double. - * - */ - void set_subsampling(double constant, double power){rate_constant = constant; rate_power = power;} + /** \brief Specifies whether the program should display information or not. + * + * @param[in] verb boolean (true = display info, false = do not display info). + * + */ + void set_verbose(bool verb = false) { verbose = verb; } public: - /** \brief Sets the mask, which is a threshold integer such that nodes in the complex that contain a number of data points which is less than or equal to - * this threshold are not displayed. - * - * @param[in] nodemask integer. - * - */ - void set_mask(int nodemask){mask = nodemask;} + /** \brief Sets the constants used to subsample the data set. These constants are + * explained in \cite Carriere17c. + * + * @param[in] constant double. + * @param[in] power double. + * + */ + void set_subsampling(double constant, double power) { + rate_constant = constant; + rate_power = power; + } public: - /** \brief Reads and stores the input point cloud. - * - * @param[in] off_file_name name of the input .OFF or .nOFF file. - * - */ - bool read_point_cloud(const std::string & off_file_name){ - - point_cloud_name = off_file_name; std::ifstream input(off_file_name); std::string line; - - char comment = '#'; - while(comment == '#'){getline(input, line); if(!line.empty() && !std::all_of(line.begin(),line.end(),isspace)) comment = line[line.find_first_not_of(' ')];} - if(std::strcmp((char*) line.c_str(),"nOFF")==0){ - comment = '#'; while(comment == '#'){getline(input, line); if(!line.empty() && !std::all_of(line.begin(),line.end(),isspace)) comment = line[line.find_first_not_of(' ')];} - std::stringstream stream(line); stream >> data_dimension; - } - else data_dimension = 3; - - comment = '#'; int numedges, numfaces, i, num; - while(comment == '#'){getline(input, line); if(!line.empty() && !std::all_of(line.begin(),line.end(),isspace)) comment = line[line.find_first_not_of(' ')];} - std::stringstream stream(line); stream >> n; stream >> numfaces; stream >> numedges; - - i = 0; while(i < n){ - getline(input, line); - if(!line.empty() && line[line.find_first_not_of(' ')] != '#' && !std::all_of(line.begin(),line.end(),isspace)){ - std::vector point; std::istringstream iss(line); - point.assign(std::istream_iterator(iss), std::istream_iterator()); - point_cloud.emplace_back(point.begin(),point.begin()+data_dimension); i++; - } - } - - i = 0; while(i < numfaces){ - getline(input, line); - if(!line.empty() && line[line.find_first_not_of(' ')] != '#' && !std::all_of(line.begin(),line.end(),isspace)){ - std::vector simplex; std::istringstream iss(line); - simplex.assign(std::istream_iterator(iss), std::istream_iterator()); - num = simplex[0]; std::vector edge(2); - for(int j = 1; j <= num; j++){ for(int k = j+1; k <= num; k++){ edge[0] = simplex[j]; edge[1] = simplex[k]; one_skeleton.push_back(edge); } } - i++; - } - } - - return input.is_open(); - } - - // ******************************************************************************************************************* - // Graphs. - // ******************************************************************************************************************* - - public: // Set graph from file. - /** \brief Creates a graph G from a file containing the edges. - * - * @param[in] graph_file_name name of the input graph file. - * The graph file contains one edge per line, - * each edge being represented by the IDs of its two nodes. - * - */ - void set_graph_from_file(const std::string & graph_file_name){ - int neighb; std::ifstream input(graph_file_name); - std::string line; int edge[2]; int n = 0; - while(std::getline(input,line)){ - std::stringstream stream(line); stream >> edge[0]; - while(stream >> neighb){edge[1] = neighb; st.insert_simplex_and_subfaces(edge);} - n++; - } - - fill_adjacency_matrix_from_st(); - - } - - - public: // Set graph from OFF file. - /** \brief Creates a graph G from the triangulation given by the input .OFF file. - * - */ - void set_graph_from_OFF(){ - - int num_edges = one_skeleton.size(); - if(num_edges > 0){ - for(int i = 0; i < num_edges; i++) st.insert_simplex_and_subfaces(one_skeleton[i]); - fill_adjacency_matrix_from_st(); - } - else std::cout << "No triangulation read in OFF file!" << std::endl; - - } - - public: // Set graph from Rips complex. - /** \brief Creates a graph G from a Rips complex. - * - * @param[in] threshold threshold value for the Rips complex. - * @param[in] distance distance used to compute the Rips complex. - * - */ - template void set_graph_from_rips(double threshold, Distance distance){ - - Rips_complex rips_complex_from_points(point_cloud, threshold, distance); - rips_complex_from_points.create_complex(st, 1); - fill_adjacency_matrix_from_st(); - - } - - - - public: // Pairwise distances. - /** \private \brief Computes all pairwise distances. - */ - template void compute_pairwise_distances(Distance ref_distance){ - - double d; std::vector zeros(n); for(int i = 0; i < n; i++) distances.push_back(zeros); - std::string distance = point_cloud_name; distance.append("_dist"); - std::ifstream input(distance.c_str(), std::ios::out | std::ios::binary); - - if(input.good()){ - if(verbose) std::cout << "Reading distances..." << std::endl; - for(int i = 0; i < n; i++){ - for (int j = i; j < n; j++){ - input.read((char*) &d,8); distances[i][j] = d; distances[j][i] = d; - } - } - input.close(); - } - - else{ - if(verbose) std::cout << "Computing distances..." << std::endl; - input.close(); std::ofstream output(distance, std::ios::out | std::ios::binary); - for(int i = 0; i < n; i++){ - int state = (int) floor( 100*(i*1.0+1)/n ) %10; - if( state == 0 && verbose) std::cout << "\r" << state << "%" << std::flush; - for (int j = i; j < n; j++){ - double dis = ref_distance(point_cloud[i],point_cloud[j]); - distances[i][j] = dis; distances[j][i] = dis; - output.write((char*) &dis,8); - } - } - output.close(); if(verbose) std::cout << std::endl; - } - - } - - public: // Automatic tuning of Rips complex. - /** \brief Creates a graph G from a Rips complex whose threshold value is automatically tuned with subsampling---see \cite Carriere17c. - * - * @param[in] distance distance between data points. - * @param[in] N number of subsampling iteration (the default reasonable value is 100, but there is no guarantee on how to choose it). - * @result delta threshold used for computing the Rips complex. - * - */ - template double set_graph_from_automatic_rips(Distance distance, int N = 100){ - - int m = floor(n/ std::exp((1+rate_power)*std::log(std::log(n)/std::log(rate_constant))) ); - m = std::min(m,n-1); - std::vector samples(m); double delta = 0; - - if(verbose) std::cout << n << " points in R^" << data_dimension << std::endl; - if(verbose) std::cout << "Subsampling " << m << " points" << std::endl; - - if(distances.size() == 0) compute_pairwise_distances(distance); - - //#pragma omp parallel for - for (int i = 0; i < N; i++){ - - SampleWithoutReplacement(n,m,samples); - double hausdorff_dist = 0; - for (int j = 0; j < n; j++){ - double mj = distances[j][samples[0]]; for (int k = 1; k < m; k++) mj = std::min(mj, distances[j][samples[k]]); - hausdorff_dist = std::max(hausdorff_dist, mj); - } - delta += hausdorff_dist/N; - - } - - if(verbose) std::cout << "delta = " << delta << std::endl; - Rips_complex rips_complex_from_points(point_cloud, delta, distance); - rips_complex_from_points.create_complex(st, 1); - fill_adjacency_matrix_from_st(); - - return delta; - - } - - - // ******************************************************************************************************************* - // Functions. - // ******************************************************************************************************************* - - public: // Set function from file. - /** \brief Creates the function f from a file containing the function values. - * - * @param[in] func_file_name name of the input function file. - * - */ - void set_function_from_file(const std::string & func_file_name){ - int vertex_id = 0; std::ifstream input(func_file_name); std::string line; double f; - while(std::getline(input,line)){ - std::stringstream stream(line); stream >> f; - func.emplace(vertex_id, f); vertex_id++; - } - functional_cover = true; - cover_name = func_file_name; - } - - public: // Set function from kth coordinate - /** \brief Creates the function f from the k-th coordinate of the point cloud P. - * - * @param[in] k coordinate to use (start at 0). - * - */ - void set_function_from_coordinate(int k){ - for(int i = 0; i < n; i++) func.emplace(i,point_cloud[i][k]); - char coordinate[100]; sprintf(coordinate, "coordinate %d", k); - functional_cover = true; - cover_name = coordinate; - } - - public: // Set function from vector. - /** \brief Creates the function f from a vector stored in memory. - * - * @param[in] function input vector of values. - * - */ - template void set_function_from_range(InputRange const & function){ - functional_cover = true; - int index = 0; for(auto v : function){func.emplace(index, v); index++;} - } - - // ******************************************************************************************************************* - // Covers. - // ******************************************************************************************************************* - - public: // Automatic tuning of resolution. - /** \brief Computes the optimal length of intervals - * (i.e. the smallest interval length avoiding discretization artifacts---see \cite Carriere17c) for a functional cover. - * - * @result reso interval length used to compute the cover. - * - */ - double set_automatic_resolution(){ - - if(!functional_cover){std::cout << "Cover needs to come from the preimages of a function." << std::endl; return 0;} - if(type != "Nerve" && type != "GIC"){std::cout << "Type of complex needs to be specified." << std::endl; return 0;} - - double reso = 0; - - if(type == "GIC"){ - for (auto simplex : st.complex_simplex_range()) { - if(st.dimension(simplex) == 1){ - std::vector vertices; - for(auto vertex : st.simplex_vertex_range(simplex)) vertices.push_back(vertex); - reso = std::max(reso, std::abs(func[vertices[0]] - func[vertices[1]])); - } - } - if(verbose) std::cout << "resolution = " << reso << std::endl; - resolution_double = reso; - } - - if(type == "Nerve"){ - for (auto simplex : st.complex_simplex_range()) { - if(st.dimension(simplex) == 1){ - std::vector vertices; - for(auto vertex : st.simplex_vertex_range(simplex)) vertices.push_back(vertex); - reso = std::max(reso, (std::abs(func[vertices[0]] - func[vertices[1]]))/gain); - } - } - if(verbose) std::cout << "resolution = " << reso << std::endl; - resolution_double = reso; - } - - return reso; - - } + /** \brief Sets the mask, which is a threshold integer such that nodes in the complex that contain a number of data + * points which is less than or equal to + * this threshold are not displayed. + * + * @param[in] nodemask integer. + * + */ + void set_mask(int nodemask) { mask = nodemask; } public: - /** \brief Sets a length of intervals from a value stored in memory. - * - * @param[in] reso length of intervals. - * - */ - void set_resolution_with_interval_length(double reso){resolution_double = reso;} - /** \brief Sets a number of intervals from a value stored in memory. - * - * @param[in] reso number of intervals. - * - */ - void set_resolution_with_interval_number(int reso){resolution_int = reso;} - /** \brief Sets a gain from a value stored in memory (default value 0.3). - * - * @param[in] g gain. - * - */ - void set_gain(double g = 0.3){gain = g;} - - - public: // Set cover with preimages of function. - /** \brief Creates a cover C from the preimages of the function f. - * - */ - void set_cover_from_function(){ - - if(resolution_double == -1 && resolution_int == -1){ std::cout << "Number and/or length of intervals not specified" << std::endl; return; } - if(gain == -1){ std::cout << "Gain not specified" << std::endl; return; } - - // Read function values and compute min and max - std::map::iterator it; - double maxf, minf; minf = std::numeric_limits::max(); maxf = std::numeric_limits::min(); - for(it = func.begin(); it != func.end(); it++){minf = std::min(minf, it->second); maxf = std::max(maxf, it->second);} - int n = func.size(); if(verbose) std::cout << "Min function value = " << minf << " and Max function value = " << maxf << std::endl; - - // Compute cover of im(f) - std::vector > intervals; int res; - - if(resolution_double == -1){ // Case we use an integer for the number of intervals. - double incr = (maxf-minf)/resolution_int; double x = minf; double alpha = (incr*gain)/(2-2*gain); - double y = minf + incr + alpha; std::pair interm(x,y); intervals.push_back(interm); - for(int i = 1; i < resolution_int-1; i++){ - x = minf + i*incr - alpha; - y = minf + (i+1)*incr + alpha; - std::pair inter(x,y); intervals.push_back(inter); - } - x = minf + (resolution_int-1)*incr - alpha; y = maxf; - std::pair interM(x,y); intervals.push_back(interM); res = intervals.size(); - if(verbose) - for(int i = 0; i < res; i++) std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" << std::endl; - } - - else{ - if(resolution_int == -1){ // Case we use a double for the length of the intervals. - double x = minf; double y = x + resolution_double; - while(y <= maxf && maxf - (y-gain*resolution_double) >= resolution_double){ - std::pair inter(x,y); intervals.push_back(inter); - x = y - gain*resolution_double; - y = x + resolution_double; - } - std::pair interM(x,maxf); intervals.push_back(interM); res = intervals.size(); - if(verbose) - for(int i = 0; i < res; i++) std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" << std::endl; - } - - else{ // Case we use an integer and a double for the length of the intervals. - double x = minf; double y = x + resolution_double; int count = 0; - while(count < resolution_int && y <= maxf && maxf - (y-gain*resolution_double) >= resolution_double){ - std::pair inter(x,y); intervals.push_back(inter); count++; - x = y - gain*resolution_double; - y = x + resolution_double; - } - res = intervals.size(); - if(verbose) - for(int i = 0; i < res; i++) std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" << std::endl; - } - } - - // Sort points according to function values - std::vector points(n); for(int i = 0; i < n; i++) points[i] = i; - std::sort(points.begin(),points.end(),Less(this->func)); - int id = 0; int pos = 0; - - for(int i = 0; i < res; i++){ - - // Find points in the preimage - std::map > prop; - std::pair inter1 = intervals[i]; - int tmp = pos; - - if(i != res-1){ - - if(i != 0){ - std::pair inter3 = intervals[i-1]; - while(func[points[tmp]] < inter3.second && tmp != n){ - prop[points[tmp]] = adjacency_matrix[points[tmp]]; - tmp++; - } - } - - std::pair inter2 = intervals[i+1]; - while(func[points[tmp]] < inter2.first && tmp != n){ - prop[points[tmp]] = adjacency_matrix[points[tmp]]; - tmp++; - } - - pos = tmp; - while(func[points[tmp]] < inter1.second && tmp != n){ - prop[points[tmp]] = adjacency_matrix[points[tmp]]; - tmp++; - } - - } - - else{ - - std::pair inter3 = intervals[i-1]; - while(func[points[tmp]] < inter3.second && tmp != n){ - prop[points[tmp]] = adjacency_matrix[points[tmp]]; - tmp++; - } - - while(tmp != n){ - prop[points[tmp]] = adjacency_matrix[points[tmp]]; - tmp++; - } - - } - - // Compute the connected components with DFS - std::map visit; if(verbose) std::cout << "Preimage of interval " << i << std::endl; - for(std::map >::iterator it = prop.begin(); it != prop.end(); it++) - visit[it->first] = false; - if (!(prop.empty())){ - for(std::map >::iterator it = prop.begin(); it != prop.end(); it++){ - if ( !(visit[it->first]) ){ - std::vector cc; cc.clear(); - dfs(prop,it->first,cc,visit); int cci = cc.size(); if(verbose) std::cout << "one CC with " << cci << " points, "; - double average_col = 0; - for(int j = 0; j < cci; j++){cover[cc[j]].push_back(id); cover_back[id].push_back(cc[j]); average_col += func_color[cc[j]]/cci;} - cover_fct[id] = i; cover_color[id] = std::pair(cci,average_col); - id++; - } - } - } - if(verbose) std::cout << std::endl; - } - - maximal_dim = id-1; - - } - - public: // Set cover from file. - /** \brief Creates the cover C from a file containing the cover elements of each point (the order has to be the same - * as in the input file!). + /** \brief Reads and stores the input point cloud. * - * @param[in] cover_file_name name of the input cover file. + * @param[in] off_file_name name of the input .OFF or .nOFF file. * */ - void set_cover_from_file(const std::string & cover_file_name){ - int vertex_id = 0; Cover_t cov; std::vector cov_elts, cov_number; - std::ifstream input(cover_file_name); std::string line; - while(std::getline(input,line)){ - cov_elts.clear(); std::stringstream stream(line); - while(stream >> cov){ - cov_elts.push_back(cov); cov_number.push_back(cov); - cover_fct[cov] = cov; cover_color[cov].second += func_color[vertex_id]; cover_color[cov].first++; cover_back[cov].push_back(vertex_id); - } - cover[vertex_id] = cov_elts; vertex_id++; - } - std::vector::iterator it; - std::sort(cov_number.begin(),cov_number.end()); it = std::unique(cov_number.begin(),cov_number.end()); - cov_number.resize(std::distance(cov_number.begin(),it)); maximal_dim = cov_number.size()-1; - for(int i = 0; i <= maximal_dim; i++) cover_color[i].second /= cover_color[i].first; - cover_name = cover_file_name; - } - - public: // Set cover from Voronoi - /** \brief Creates the cover C from the Voronoï cells of a subsampling of the point cloud. + bool read_point_cloud(const std::string& off_file_name) { + point_cloud_name = off_file_name; + std::ifstream input(off_file_name); + std::string line; + + char comment = '#'; + while (comment == '#') { + getline(input, line); + if (!line.empty() && !std::all_of(line.begin(), line.end(), isspace)) comment = line[line.find_first_not_of(' ')]; + } + if (std::strcmp((char*)line.c_str(), "nOFF") == 0) { + comment = '#'; + while (comment == '#') { + getline(input, line); + if (!line.empty() && !std::all_of(line.begin(), line.end(), isspace)) + comment = line[line.find_first_not_of(' ')]; + } + std::stringstream stream(line); + stream >> data_dimension; + } else + data_dimension = 3; + + comment = '#'; + int numedges, numfaces, i, num; + while (comment == '#') { + getline(input, line); + if (!line.empty() && !std::all_of(line.begin(), line.end(), isspace)) comment = line[line.find_first_not_of(' ')]; + } + std::stringstream stream(line); + stream >> n; + stream >> numfaces; + stream >> numedges; + + i = 0; + while (i < n) { + getline(input, line); + if (!line.empty() && line[line.find_first_not_of(' ')] != '#' && + !std::all_of(line.begin(), line.end(), isspace)) { + std::vector point; + std::istringstream iss(line); + point.assign(std::istream_iterator(iss), std::istream_iterator()); + point_cloud.emplace_back(point.begin(), point.begin() + data_dimension); + i++; + } + } + + i = 0; + while (i < numfaces) { + getline(input, line); + if (!line.empty() && line[line.find_first_not_of(' ')] != '#' && + !std::all_of(line.begin(), line.end(), isspace)) { + std::vector simplex; + std::istringstream iss(line); + simplex.assign(std::istream_iterator(iss), std::istream_iterator()); + num = simplex[0]; + std::vector edge(2); + for (int j = 1; j <= num; j++) { + for (int k = j + 1; k <= num; k++) { + edge[0] = simplex[j]; + edge[1] = simplex[k]; + one_skeleton.push_back(edge); + } + } + i++; + } + } + + return input.is_open(); + } + + // ******************************************************************************************************************* + // Graphs. + // ******************************************************************************************************************* + + public: // Set graph from file. + /** \brief Creates a graph G from a file containing the edges. + * + * @param[in] graph_file_name name of the input graph file. + * The graph file contains one edge per line, + * each edge being represented by the IDs of its two nodes. + * + */ + void set_graph_from_file(const std::string& graph_file_name) { + int neighb; + std::ifstream input(graph_file_name); + std::string line; + int edge[2]; + int n = 0; + while (std::getline(input, line)) { + std::stringstream stream(line); + stream >> edge[0]; + while (stream >> neighb) { + edge[1] = neighb; + st.insert_simplex_and_subfaces(edge); + } + n++; + } + + fill_adjacency_matrix_from_st(); + } + + public: // Set graph from OFF file. + /** \brief Creates a graph G from the triangulation given by the input .OFF file. + * + */ + void set_graph_from_OFF() { + int num_edges = one_skeleton.size(); + if (num_edges > 0) { + for (int i = 0; i < num_edges; i++) st.insert_simplex_and_subfaces(one_skeleton[i]); + fill_adjacency_matrix_from_st(); + } else + std::cout << "No triangulation read in OFF file!" << std::endl; + } + + public: // Set graph from Rips complex. + /** \brief Creates a graph G from a Rips complex. + * + * @param[in] threshold threshold value for the Rips complex. + * @param[in] distance distance used to compute the Rips complex. + * + */ + template + void set_graph_from_rips(double threshold, Distance distance) { + Rips_complex rips_complex_from_points(point_cloud, threshold, distance); + rips_complex_from_points.create_complex(st, 1); + fill_adjacency_matrix_from_st(); + } + + public: // Pairwise distances. + /** \private \brief Computes all pairwise distances. + */ + template + void compute_pairwise_distances(Distance ref_distance) { + double d; + std::vector zeros(n); + for (int i = 0; i < n; i++) distances.push_back(zeros); + std::string distance = point_cloud_name; + distance.append("_dist"); + std::ifstream input(distance.c_str(), std::ios::out | std::ios::binary); + + if (input.good()) { + if (verbose) std::cout << "Reading distances..." << std::endl; + for (int i = 0; i < n; i++) { + for (int j = i; j < n; j++) { + input.read((char*)&d, 8); + distances[i][j] = d; + distances[j][i] = d; + } + } + input.close(); + } + + else { + if (verbose) std::cout << "Computing distances..." << std::endl; + input.close(); + std::ofstream output(distance, std::ios::out | std::ios::binary); + for (int i = 0; i < n; i++) { + int state = (int)floor(100 * (i * 1.0 + 1) / n) % 10; + if (state == 0 && verbose) std::cout << "\r" << state << "%" << std::flush; + for (int j = i; j < n; j++) { + double dis = ref_distance(point_cloud[i], point_cloud[j]); + distances[i][j] = dis; + distances[j][i] = dis; + output.write((char*)&dis, 8); + } + } + output.close(); + if (verbose) std::cout << std::endl; + } + } + + public: // Automatic tuning of Rips complex. + /** \brief Creates a graph G from a Rips complex whose threshold value is automatically tuned with subsampling---see + * \cite Carriere17c. * - * @param[in] distance distance between the points. - * @param[in] m number of points in the subsample. + * @param[in] distance distance between data points. + * @param[in] N number of subsampling iteration (the default reasonable value is 100, but there is no guarantee on + * how to choose it). + * @result delta threshold used for computing the Rips complex. * */ - template void set_cover_from_Voronoi(Distance distance, int m = 100){ - - voronoi_subsamples.resize(m); SampleWithoutReplacement(n,m,voronoi_subsamples); - if(distances.size() == 0) compute_pairwise_distances(distance); - std::vector mindist(n); for(int j = 0; j < n; j++) mindist[j] = std::numeric_limits::max(); - - // Compute the geodesic distances to subsamples with Dijkstra - for(int i = 0; i < m; i++){ - if(verbose) std::cout << "Computing geodesic distances to seed " << i << "..." << std::endl; - int seed = voronoi_subsamples[i]; - std::vector dist(n); std::vector process(n); - for(int j = 0; j < n; j++){ dist[j] = std::numeric_limits::max(); process[j] = j; } - dist[seed] = 0; int curr_size = process.size(); int min_point, min_index; double min_dist; - std::vector neighbors; int num_neighbors; - - while(curr_size > 0){ - min_dist = std::numeric_limits::max(); min_index = -1; min_point = -1; - for(int j = 0; j < curr_size; j++){ - if(dist[process[j]] < min_dist){ - min_point = process[j]; min_dist = dist[process[j]]; min_index = j; - } - } - assert(min_index != -1); process.erase(process.begin() + min_index); - assert(min_point != -1); neighbors = adjacency_matrix[min_point]; num_neighbors = neighbors.size(); - for(int j = 0; j < num_neighbors; j++){ - double d = dist[min_point] + distances[min_point][neighbors[j]]; - dist[neighbors[j]] = std::min(dist[neighbors[j]], d); - } - curr_size = process.size(); - } - - for(int j = 0; j < n; j++) - if(mindist[j] > dist[j]){ - mindist[j] = dist[j]; - if(cover[j].size() == 0) cover[j].push_back(i); - else cover[j][0] = i; - } - } - - for(int i = 0; i < n; i++){ cover_back[cover[i][0]].push_back(i); cover_color[cover[i][0]].second += func_color[i]; cover_color[cover[i][0]].first++; } - for(int i = 0; i < m; i++) cover_color[i].second /= cover_color[i].first; - maximal_dim = m-1; cover_name = "Voronoi"; - - } - - public: // return subset of data corresponding to a node - /** \brief Returns the data subset corresponding to a specific node of the created complex. + template + double set_graph_from_automatic_rips(Distance distance, int N = 100) { + int m = floor(n / std::exp((1 + rate_power) * std::log(std::log(n) / std::log(rate_constant)))); + m = std::min(m, n - 1); + std::vector samples(m); + double delta = 0; + + if (verbose) std::cout << n << " points in R^" << data_dimension << std::endl; + if (verbose) std::cout << "Subsampling " << m << " points" << std::endl; + + if (distances.size() == 0) compute_pairwise_distances(distance); + + //#pragma omp parallel for + for (int i = 0; i < N; i++) { + SampleWithoutReplacement(n, m, samples); + double hausdorff_dist = 0; + for (int j = 0; j < n; j++) { + double mj = distances[j][samples[0]]; + for (int k = 1; k < m; k++) mj = std::min(mj, distances[j][samples[k]]); + hausdorff_dist = std::max(hausdorff_dist, mj); + } + delta += hausdorff_dist / N; + } + + if (verbose) std::cout << "delta = " << delta << std::endl; + Rips_complex rips_complex_from_points(point_cloud, delta, distance); + rips_complex_from_points.create_complex(st, 1); + fill_adjacency_matrix_from_st(); + + return delta; + } + + // ******************************************************************************************************************* + // Functions. + // ******************************************************************************************************************* + + public: // Set function from file. + /** \brief Creates the function f from a file containing the function values. + * + * @param[in] func_file_name name of the input function file. + * + */ + void set_function_from_file(const std::string& func_file_name) { + int vertex_id = 0; + std::ifstream input(func_file_name); + std::string line; + double f; + while (std::getline(input, line)) { + std::stringstream stream(line); + stream >> f; + func.emplace(vertex_id, f); + vertex_id++; + } + functional_cover = true; + cover_name = func_file_name; + } + + public: // Set function from kth coordinate + /** \brief Creates the function f from the k-th coordinate of the point cloud P. + * + * @param[in] k coordinate to use (start at 0). + * + */ + void set_function_from_coordinate(int k) { + for (int i = 0; i < n; i++) func.emplace(i, point_cloud[i][k]); + char coordinate[100]; + sprintf(coordinate, "coordinate %d", k); + functional_cover = true; + cover_name = coordinate; + } + + public: // Set function from vector. + /** \brief Creates the function f from a vector stored in memory. + * + * @param[in] function input vector of values. + * + */ + template + void set_function_from_range(InputRange const& function) { + functional_cover = true; + int index = 0; + for (auto v : function) { + func.emplace(index, v); + index++; + } + } + + // ******************************************************************************************************************* + // Covers. + // ******************************************************************************************************************* + + public: // Automatic tuning of resolution. + /** \brief Computes the optimal length of intervals + * (i.e. the smallest interval length avoiding discretization artifacts---see \cite Carriere17c) for a functional + * cover. + * + * @result reso interval length used to compute the cover. + * + */ + double set_automatic_resolution() { + if (!functional_cover) { + std::cout << "Cover needs to come from the preimages of a function." << std::endl; + return 0; + } + if (type != "Nerve" && type != "GIC") { + std::cout << "Type of complex needs to be specified." << std::endl; + return 0; + } + + double reso = 0; + + if (type == "GIC") { + for (auto simplex : st.complex_simplex_range()) { + if (st.dimension(simplex) == 1) { + std::vector vertices; + for (auto vertex : st.simplex_vertex_range(simplex)) vertices.push_back(vertex); + reso = std::max(reso, std::abs(func[vertices[0]] - func[vertices[1]])); + } + } + if (verbose) std::cout << "resolution = " << reso << std::endl; + resolution_double = reso; + } + + if (type == "Nerve") { + for (auto simplex : st.complex_simplex_range()) { + if (st.dimension(simplex) == 1) { + std::vector vertices; + for (auto vertex : st.simplex_vertex_range(simplex)) vertices.push_back(vertex); + reso = std::max(reso, (std::abs(func[vertices[0]] - func[vertices[1]])) / gain); + } + } + if (verbose) std::cout << "resolution = " << reso << std::endl; + resolution_double = reso; + } + + return reso; + } + + public: + /** \brief Sets a length of intervals from a value stored in memory. * - * @param[in] c ID of the node. - * @result cover_back(c) vector of IDs of data points. + * @param[in] reso length of intervals. * */ - const std::vector & subpopulation(Cover_t c){ - return cover_back[c]; - } - - // ******************************************************************************************************************* - // Visualization. - // ******************************************************************************************************************* - - public: // Set color from file. - /** \brief Computes the function used to color the nodes of the simplicial complex from a file containing the function values. - * - * @param[in] color_file_name name of the input color 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.emplace(vertex_id, f); vertex_id++; - } - color_name = color_file_name; - } - - public: // Set color from kth coordinate - /** \brief Computes the function used to color the nodes of the simplicial complex from the k-th coordinate. - * - * @param[in] k coordinate to use (start at 0). - * - */ - void set_color_from_coordinate(int k = 0){ - for(int i = 0; i < n; i++) func_color.emplace(i, point_cloud[i][k]); - color_name = "coordinate "; color_name.append(std::to_string(k)); - } - - public: // Set color from vector. - /** \brief Computes the function used to color the nodes of the simplicial complex from a vector stored in memory. - * - * @param[in] color input vector of values. - * - */ - void set_color_from_vector(std::vector color){ - for(unsigned int i = 0; i < color.size(); i++) func_color.emplace(i, color[i]); - } - - public: // Create a .dot file that can be compiled with neato to produce a .pdf file. - /** \brief Creates a .dot file called SC.dot for neato (part of the graphviz package) once the simplicial complex is computed to get a visualization - * of its 1-skeleton in a .pdf file. - */ - void plot_DOT(){ - char mapp[11] = "SC.dot"; std::ofstream graphic(mapp); graphic << "graph GIC {" << std::endl; - double maxv, minv; maxv = std::numeric_limits::min(); minv = std::numeric_limits::max(); - for (std::map >::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 nodes; nodes.clear(); - for (std::map >::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(); - std::cout << "SC.dot generated. It can be visualized with e.g. neato." << std::endl; - } - - public: // Create a .txt file that can be compiled with KeplerMapper. - /** \brief Creates a .txt file called SC.txt describing the 1-skeleton, which can then be plotted with e.g. KeplerMapper. - */ - void write_info(){ - - int num_simplices = simplices.size(); int num_edges = 0; - char mapp[11] = "SC.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 << point_cloud_name << std::endl; - graphic << cover_name << std::endl; - graphic << color_name << std::endl; - graphic << resolution_double << " " << gain << std::endl; - graphic << cover_color.size() << " " << num_edges << std::endl; - - for (std::map >::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(); - std::cout << "SC.txt generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox." << std::endl; - } - - - public: // Create a .off file that can be visualized (e.g. with Geomview). - /** \brief Creates a .off file called SC.off for 3D visualization, which contains the 2-skeleton of the GIC. - * This function assumes that the cover has been computed with Voronoi. If data points are in 1D or 2D, - * the remaining coordinates of the points embedded in 3D are set to 0. - */ - void plot_OFF(){ - - assert(cover_name == "Voronoi"); - char gic[11] = "SC.off"; std::ofstream graphic(gic); - graphic << "OFF" << std::endl; int m = voronoi_subsamples.size(); int numedges = 0; int numfaces = 0; - std::vector > edges, faces; int numsimplices = simplices.size(); - for (int i = 0; i < numsimplices; i++) { - if(simplices[i].size() == 2){ numedges++; edges.push_back(simplices[i]); } - if(simplices[i].size() == 3){ numfaces++; faces.push_back(simplices[i]); } - } - graphic << m << " " << numedges + numfaces << std::endl; - for(int i = 0; i < m; i++){ - if(data_dimension <= 3){ - for(int j = 0; j < data_dimension; j++) graphic << point_cloud[voronoi_subsamples[i]][j] << " "; - for(int j = data_dimension; j < 3; j++) graphic << 0 << " "; - graphic << std::endl; - } - else for(int j = 0; j < 3; j++) graphic << point_cloud[voronoi_subsamples[i]][j] << " "; - } - for(int i = 0; i < numedges; i++) graphic << 2 << " " << edges[i][0] << " " << edges[i][1] << std::endl; - for(int i = 0; i < numfaces; i++) graphic << 3 << " " << faces[i][0] << " " << faces[i][1] << " " << faces[i][2] << std::endl; - graphic.close(); - std::cout << "SC.off generated. It can be visualized with e.g. geomview." << std::endl; - - } - - // ******************************************************************************************************************* - // ******************************************************************************************************************* - + void set_resolution_with_interval_length(double reso) { resolution_double = reso; } + /** \brief Sets a number of intervals from a value stored in memory. + * + * @param[in] reso number of intervals. + * + */ + void set_resolution_with_interval_number(int reso) { resolution_int = reso; } + /** \brief Sets a gain from a value stored in memory (default value 0.3). + * + * @param[in] g gain. + * + */ + void set_gain(double g = 0.3) { gain = g; } + + public: // Set cover with preimages of function. + /** \brief Creates a cover C from the preimages of the function f. + * + */ + void set_cover_from_function() { + if (resolution_double == -1 && resolution_int == -1) { + std::cout << "Number and/or length of intervals not specified" << std::endl; + return; + } + if (gain == -1) { + std::cout << "Gain not specified" << std::endl; + return; + } + + // Read function values and compute min and max + std::map::iterator it; + double maxf, minf; + minf = std::numeric_limits::max(); + maxf = std::numeric_limits::min(); + for (it = func.begin(); it != func.end(); it++) { + minf = std::min(minf, it->second); + maxf = std::max(maxf, it->second); + } + int n = func.size(); + if (verbose) std::cout << "Min function value = " << minf << " and Max function value = " << maxf << std::endl; + + // Compute cover of im(f) + std::vector > intervals; + int res; + + if (resolution_double == -1) { // Case we use an integer for the number of intervals. + double incr = (maxf - minf) / resolution_int; + double x = minf; + double alpha = (incr * gain) / (2 - 2 * gain); + double y = minf + incr + alpha; + std::pair interm(x, y); + intervals.push_back(interm); + for (int i = 1; i < resolution_int - 1; i++) { + x = minf + i * incr - alpha; + y = minf + (i + 1) * incr + alpha; + std::pair inter(x, y); + intervals.push_back(inter); + } + x = minf + (resolution_int - 1) * incr - alpha; + y = maxf; + std::pair interM(x, y); + intervals.push_back(interM); + res = intervals.size(); + if (verbose) + for (int i = 0; i < res; i++) + std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" + << std::endl; + } + + else { + if (resolution_int == -1) { // Case we use a double for the length of the intervals. + double x = minf; + double y = x + resolution_double; + while (y <= maxf && maxf - (y - gain * resolution_double) >= resolution_double) { + std::pair inter(x, y); + intervals.push_back(inter); + x = y - gain * resolution_double; + y = x + resolution_double; + } + std::pair interM(x, maxf); + intervals.push_back(interM); + res = intervals.size(); + if (verbose) + for (int i = 0; i < res; i++) + std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" + << std::endl; + } + + else { // Case we use an integer and a double for the length of the intervals. + double x = minf; + double y = x + resolution_double; + int count = 0; + while (count < resolution_int && y <= maxf && maxf - (y - gain * resolution_double) >= resolution_double) { + std::pair inter(x, y); + intervals.push_back(inter); + count++; + x = y - gain * resolution_double; + y = x + resolution_double; + } + res = intervals.size(); + if (verbose) + for (int i = 0; i < res; i++) + std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" + << std::endl; + } + } + + // Sort points according to function values + std::vector points(n); + for (int i = 0; i < n; i++) points[i] = i; + std::sort(points.begin(), points.end(), Less(this->func)); + int id = 0; + int pos = 0; + + for (int i = 0; i < res; i++) { + // Find points in the preimage + std::map > prop; + std::pair inter1 = intervals[i]; + int tmp = pos; + + if (i != res - 1) { + if (i != 0) { + std::pair inter3 = intervals[i - 1]; + while (func[points[tmp]] < inter3.second && tmp != n) { + prop[points[tmp]] = adjacency_matrix[points[tmp]]; + tmp++; + } + } + + std::pair inter2 = intervals[i + 1]; + while (func[points[tmp]] < inter2.first && tmp != n) { + prop[points[tmp]] = adjacency_matrix[points[tmp]]; + tmp++; + } + + pos = tmp; + while (func[points[tmp]] < inter1.second && tmp != n) { + prop[points[tmp]] = adjacency_matrix[points[tmp]]; + tmp++; + } + + } + + else { + std::pair inter3 = intervals[i - 1]; + while (func[points[tmp]] < inter3.second && tmp != n) { + prop[points[tmp]] = adjacency_matrix[points[tmp]]; + tmp++; + } + + while (tmp != n) { + prop[points[tmp]] = adjacency_matrix[points[tmp]]; + tmp++; + } + } + + // Compute the connected components with DFS + std::map visit; + if (verbose) std::cout << "Preimage of interval " << i << std::endl; + for (std::map >::iterator it = prop.begin(); it != prop.end(); it++) + visit[it->first] = false; + if (!(prop.empty())) { + for (std::map >::iterator it = prop.begin(); it != prop.end(); it++) { + if (!(visit[it->first])) { + std::vector cc; + cc.clear(); + dfs(prop, it->first, cc, visit); + int cci = cc.size(); + if (verbose) std::cout << "one CC with " << cci << " points, "; + double average_col = 0; + for (int j = 0; j < cci; j++) { + cover[cc[j]].push_back(id); + cover_back[id].push_back(cc[j]); + average_col += func_color[cc[j]] / cci; + } + cover_fct[id] = i; + cover_color[id] = std::pair(cci, average_col); + id++; + } + } + } + if (verbose) std::cout << std::endl; + } + + maximal_dim = id - 1; + } + + public: // Set cover from file. + /** \brief Creates the cover C from a file containing the cover elements of each point (the order has to be the same + * as in the input file!). + * + * @param[in] cover_file_name name of the input cover file. + * + */ + void set_cover_from_file(const std::string& cover_file_name) { + int vertex_id = 0; + Cover_t cov; + std::vector cov_elts, cov_number; + std::ifstream input(cover_file_name); + std::string line; + while (std::getline(input, line)) { + cov_elts.clear(); + std::stringstream stream(line); + while (stream >> cov) { + cov_elts.push_back(cov); + cov_number.push_back(cov); + cover_fct[cov] = cov; + cover_color[cov].second += func_color[vertex_id]; + cover_color[cov].first++; + cover_back[cov].push_back(vertex_id); + } + cover[vertex_id] = cov_elts; + vertex_id++; + } + std::vector::iterator it; + std::sort(cov_number.begin(), cov_number.end()); + it = std::unique(cov_number.begin(), cov_number.end()); + cov_number.resize(std::distance(cov_number.begin(), it)); + maximal_dim = cov_number.size() - 1; + for (int i = 0; i <= maximal_dim; i++) cover_color[i].second /= cover_color[i].first; + cover_name = cover_file_name; + } + + public: // Set cover from Voronoi + /** \brief Creates the cover C from the Voronoï cells of a subsampling of the point cloud. + * + * @param[in] distance distance between the points. + * @param[in] m number of points in the subsample. + * + */ + template + void set_cover_from_Voronoi(Distance distance, int m = 100) { + voronoi_subsamples.resize(m); + SampleWithoutReplacement(n, m, voronoi_subsamples); + if (distances.size() == 0) compute_pairwise_distances(distance); + std::vector mindist(n); + for (int j = 0; j < n; j++) mindist[j] = std::numeric_limits::max(); + + // Compute the geodesic distances to subsamples with Dijkstra + for (int i = 0; i < m; i++) { + if (verbose) std::cout << "Computing geodesic distances to seed " << i << "..." << std::endl; + int seed = voronoi_subsamples[i]; + std::vector dist(n); + std::vector process(n); + for (int j = 0; j < n; j++) { + dist[j] = std::numeric_limits::max(); + process[j] = j; + } + dist[seed] = 0; + int curr_size = process.size(); + int min_point, min_index; + double min_dist; + std::vector neighbors; + int num_neighbors; + + while (curr_size > 0) { + min_dist = std::numeric_limits::max(); + min_index = -1; + min_point = -1; + for (int j = 0; j < curr_size; j++) { + if (dist[process[j]] < min_dist) { + min_point = process[j]; + min_dist = dist[process[j]]; + min_index = j; + } + } + assert(min_index != -1); + process.erase(process.begin() + min_index); + assert(min_point != -1); + neighbors = adjacency_matrix[min_point]; + num_neighbors = neighbors.size(); + for (int j = 0; j < num_neighbors; j++) { + double d = dist[min_point] + distances[min_point][neighbors[j]]; + dist[neighbors[j]] = std::min(dist[neighbors[j]], d); + } + curr_size = process.size(); + } + + for (int j = 0; j < n; j++) + if (mindist[j] > dist[j]) { + mindist[j] = dist[j]; + if (cover[j].size() == 0) + cover[j].push_back(i); + else + cover[j][0] = i; + } + } + + for (int i = 0; i < n; i++) { + cover_back[cover[i][0]].push_back(i); + cover_color[cover[i][0]].second += func_color[i]; + cover_color[cover[i][0]].first++; + } + for (int i = 0; i < m; i++) cover_color[i].second /= cover_color[i].first; + maximal_dim = m - 1; + cover_name = "Voronoi"; + } + + public: // return subset of data corresponding to a node + /** \brief Returns the data subset corresponding to a specific node of the created complex. + * + * @param[in] c ID of the node. + * @result cover_back(c) vector of IDs of data points. + * + */ + const std::vector& subpopulation(Cover_t c) { return cover_back[c]; } + + // ******************************************************************************************************************* + // Visualization. + // ******************************************************************************************************************* + + public: // Set color from file. + /** \brief Computes the function used to color the nodes of the simplicial complex from a file containing the function + * values. + * + * @param[in] color_file_name name of the input color 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.emplace(vertex_id, f); + vertex_id++; + } + color_name = color_file_name; + } + + public: // Set color from kth coordinate + /** \brief Computes the function used to color the nodes of the simplicial complex from the k-th coordinate. + * + * @param[in] k coordinate to use (start at 0). + * + */ + void set_color_from_coordinate(int k = 0) { + for (int i = 0; i < n; i++) func_color.emplace(i, point_cloud[i][k]); + color_name = "coordinate "; + color_name.append(std::to_string(k)); + } + + public: // Set color from vector. + /** \brief Computes the function used to color the nodes of the simplicial complex from a vector stored in memory. + * + * @param[in] color input vector of values. + * + */ + void set_color_from_vector(std::vector color) { + for (unsigned int i = 0; i < color.size(); i++) func_color.emplace(i, color[i]); + } + + public: // Create a .dot file that can be compiled with neato to produce a .pdf file. + /** \brief Creates a .dot file called SC.dot for neato (part of the graphviz package) once the simplicial complex is + * computed to get a visualization + * of its 1-skeleton in a .pdf file. + */ + void plot_DOT() { + char mapp[11] = "SC.dot"; + std::ofstream graphic(mapp); + graphic << "graph GIC {" << std::endl; + double maxv, minv; + maxv = std::numeric_limits::min(); + minv = std::numeric_limits::max(); + for (std::map >::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 nodes; + nodes.clear(); + for (std::map >::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(); + std::cout << "SC.dot generated. It can be visualized with e.g. neato." << std::endl; + } + + public: // Create a .txt file that can be compiled with KeplerMapper. + /** \brief Creates a .txt file called SC.txt describing the 1-skeleton, which can then be plotted with e.g. + * KeplerMapper. + */ + void write_info() { + int num_simplices = simplices.size(); + int num_edges = 0; + char mapp[11] = "SC.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 << point_cloud_name << std::endl; + graphic << cover_name << std::endl; + graphic << color_name << std::endl; + graphic << resolution_double << " " << gain << std::endl; + graphic << cover_color.size() << " " << num_edges << std::endl; + + for (std::map >::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(); + std::cout << "SC.txt generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox." + << std::endl; + } + + public: // Create a .off file that can be visualized (e.g. with Geomview). + /** \brief Creates a .off file called SC.off for 3D visualization, which contains the 2-skeleton of the GIC. + * This function assumes that the cover has been computed with Voronoi. If data points are in 1D or 2D, + * the remaining coordinates of the points embedded in 3D are set to 0. + */ + void plot_OFF() { + assert(cover_name == "Voronoi"); + char gic[11] = "SC.off"; + std::ofstream graphic(gic); + graphic << "OFF" << std::endl; + int m = voronoi_subsamples.size(); + int numedges = 0; + int numfaces = 0; + std::vector > edges, faces; + int numsimplices = simplices.size(); + for (int i = 0; i < numsimplices; i++) { + if (simplices[i].size() == 2) { + numedges++; + edges.push_back(simplices[i]); + } + if (simplices[i].size() == 3) { + numfaces++; + faces.push_back(simplices[i]); + } + } + graphic << m << " " << numedges + numfaces << std::endl; + for (int i = 0; i < m; i++) { + if (data_dimension <= 3) { + for (int j = 0; j < data_dimension; j++) graphic << point_cloud[voronoi_subsamples[i]][j] << " "; + for (int j = data_dimension; j < 3; j++) graphic << 0 << " "; + graphic << std::endl; + } else + for (int j = 0; j < 3; j++) graphic << point_cloud[voronoi_subsamples[i]][j] << " "; + } + for (int i = 0; i < numedges; i++) graphic << 2 << " " << edges[i][0] << " " << edges[i][1] << std::endl; + for (int i = 0; i < numfaces; i++) + graphic << 3 << " " << faces[i][0] << " " << faces[i][1] << " " << faces[i][2] << std::endl; + graphic.close(); + std::cout << "SC.off generated. It can be visualized with e.g. geomview." << std::endl; + } + + // ******************************************************************************************************************* + // ******************************************************************************************************************* public: - /** \brief Creates the simplicial complex. - * - * @param[in] complex SimplicialComplexForRips to be created. - * - */ - template - void create_complex(SimplicialComplexForRips & complex) { - unsigned int dimension = 0; - for(auto const& simplex : simplices){ - complex.insert_simplex_and_subfaces(simplex); - if(dimension < simplex.size()-1) dimension = simplex.size()-1; - } - complex.set_dimension(dimension); - } + /** \brief Creates the simplicial complex. + * + * @param[in] complex SimplicialComplexForRips to be created. + * + */ + template + void create_complex(SimplicialComplexForRips& complex) { + unsigned int dimension = 0; + for (auto const& simplex : simplices) { + complex.insert_simplex_and_subfaces(simplex); + if (dimension < simplex.size() - 1) dimension = simplex.size() - 1; + } + complex.set_dimension(dimension); + } public: - /** \brief Computes the simplices of the simplicial complex. - */ - void find_simplices(){ - - if(type != "Nerve" && type != "GIC"){std::cout << "Type of complex needs to be specified." << std::endl; return;} - - if(type == "Nerve"){ - for(std::map >::iterator it = cover.begin(); it!=cover.end(); it++) simplices.push_back(it->second); - std::vector >::iterator it; - std::sort(simplices.begin(),simplices.end()); it = std::unique(simplices.begin(),simplices.end()); - simplices.resize(std::distance(simplices.begin(),it)); - } - - if(type == "GIC"){ - - if(functional_cover){ - - // Computes the simplices in the GIC by looking at all the edges of the graph and adding the - // corresponding edges in the GIC if the images of the endpoints belong to consecutive intervals. - - if (gain >= 0.5) - throw std::invalid_argument("the output of this function is correct ONLY if the cover is minimal, i.e. the gain is less than 0.5."); - - int v1, v2; - - // Loop on all points. - for(std::map >::iterator it = cover.begin(); it != cover.end(); it++){ - - int vid = it->first; std::vector neighbors = adjacency_matrix[vid]; int num_neighb = neighbors.size(); - - // 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 node(1); node[0] = v1; simplices.push_back(node); - - // Loop on neighbors. - for(int i = 0; i < num_neighb; i++){ - - int neighb = neighbors[i]; - - // Find cover of neighbor (neighb). - if(cover[neighb].size() == 2) v2 = std::max(cover[neighb][0],cover[neighb][1]); else v2 = cover[neighb][0]; - - // If neighbor is in next interval, add edge. - if(cover_fct[v2] == cover_fct[v1] + 1){ - std::vector edge(2); edge[0] = v1; edge[1] = v2; - simplices.push_back(edge); - } - } - } - std::vector >::iterator it; - std::sort(simplices.begin(),simplices.end()); it = std::unique(simplices.begin(),simplices.end()); - simplices.resize(std::distance(simplices.begin(),it)); - - } - - else{ - - // Find IDs of edges to remove - std::vector simplex_to_remove; int simplex_id = 0; - for (auto simplex : st.complex_simplex_range()) { - if(st.dimension(simplex) == 1){ - std::vector > comp; - for(auto vertex : st.simplex_vertex_range(simplex)) comp.push_back(cover[vertex]); - if(comp[0].size() == 1 && comp[0] == comp[1]) simplex_to_remove.push_back(simplex_id); - } - simplex_id++; - } - - // Remove edges - if(simplex_to_remove.size() > 1){ - int current_id = 1; - auto simplex = st.complex_simplex_range().begin(); int num_rem = 0; - for(int i = 0; i < simplex_id-1; i++){ - int j = i+1; auto simplex_tmp = simplex; simplex_tmp++; - if(j == simplex_to_remove[current_id]){st.remove_maximal_simplex(*simplex_tmp); current_id++; num_rem++;} - else simplex++; - } simplex = st.complex_simplex_range().begin(); - for(int i = 0; i < simplex_to_remove[0]; i++) simplex++; st.remove_maximal_simplex(*simplex); - } - - // Build the Simplex Tree corresponding to the graph - st.expansion(maximal_dim); - - // Find simplices of GIC - simplices.clear(); - for (auto simplex : st.complex_simplex_range()) { - if(!st.has_children(simplex)){ - - std::vector simplx; - for (auto vertex : st.simplex_vertex_range(simplex)){ - unsigned int sz = cover[vertex].size(); - for(unsigned int i = 0; i < sz; i++){ - simplx.push_back(cover[vertex][i]); - } - } - - std::sort(simplx.begin(),simplx.end()); std::vector::iterator it = std::unique(simplx.begin(),simplx.end()); - simplx.resize(std::distance(simplx.begin(),it)); - simplices.push_back(simplx); - - } - } - std::vector >::iterator it; - std::sort(simplices.begin(),simplices.end()); it = std::unique(simplices.begin(),simplices.end()); - simplices.resize(std::distance(simplices.begin(),it)); - - } - } - - } - + /** \brief Computes the simplices of the simplicial complex. + */ + void find_simplices() { + if (type != "Nerve" && type != "GIC") { + std::cout << "Type of complex needs to be specified." << std::endl; + return; + } + + if (type == "Nerve") { + for (std::map >::iterator it = cover.begin(); it != cover.end(); it++) + simplices.push_back(it->second); + std::vector >::iterator it; + std::sort(simplices.begin(), simplices.end()); + it = std::unique(simplices.begin(), simplices.end()); + simplices.resize(std::distance(simplices.begin(), it)); + } + + if (type == "GIC") { + if (functional_cover) { + // Computes the simplices in the GIC by looking at all the edges of the graph and adding the + // corresponding edges in the GIC if the images of the endpoints belong to consecutive intervals. + + if (gain >= 0.5) + throw std::invalid_argument( + "the output of this function is correct ONLY if the cover is minimal, i.e. the gain is less than 0.5."); + + int v1, v2; + + // Loop on all points. + for (std::map >::iterator it = cover.begin(); it != cover.end(); it++) { + int vid = it->first; + std::vector neighbors = adjacency_matrix[vid]; + int num_neighb = neighbors.size(); + + // 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 node(1); + node[0] = v1; + simplices.push_back(node); + + // Loop on neighbors. + for (int i = 0; i < num_neighb; i++) { + int neighb = neighbors[i]; + + // Find cover of neighbor (neighb). + if (cover[neighb].size() == 2) + v2 = std::max(cover[neighb][0], cover[neighb][1]); + else + v2 = cover[neighb][0]; + + // If neighbor is in next interval, add edge. + if (cover_fct[v2] == cover_fct[v1] + 1) { + std::vector edge(2); + edge[0] = v1; + edge[1] = v2; + simplices.push_back(edge); + } + } + } + std::vector >::iterator it; + std::sort(simplices.begin(), simplices.end()); + it = std::unique(simplices.begin(), simplices.end()); + simplices.resize(std::distance(simplices.begin(), it)); + + } + + else { + // Find IDs of edges to remove + std::vector simplex_to_remove; + int simplex_id = 0; + for (auto simplex : st.complex_simplex_range()) { + if (st.dimension(simplex) == 1) { + std::vector > comp; + for (auto vertex : st.simplex_vertex_range(simplex)) comp.push_back(cover[vertex]); + if (comp[0].size() == 1 && comp[0] == comp[1]) simplex_to_remove.push_back(simplex_id); + } + simplex_id++; + } + + // Remove edges + if (simplex_to_remove.size() > 1) { + int current_id = 1; + auto simplex = st.complex_simplex_range().begin(); + int num_rem = 0; + for (int i = 0; i < simplex_id - 1; i++) { + int j = i + 1; + auto simplex_tmp = simplex; + simplex_tmp++; + if (j == simplex_to_remove[current_id]) { + st.remove_maximal_simplex(*simplex_tmp); + current_id++; + num_rem++; + } else + simplex++; + } + simplex = st.complex_simplex_range().begin(); + for (int i = 0; i < simplex_to_remove[0]; i++) simplex++; + st.remove_maximal_simplex(*simplex); + } + + // Build the Simplex Tree corresponding to the graph + st.expansion(maximal_dim); + + // Find simplices of GIC + simplices.clear(); + for (auto simplex : st.complex_simplex_range()) { + if (!st.has_children(simplex)) { + std::vector simplx; + for (auto vertex : st.simplex_vertex_range(simplex)) { + unsigned int sz = cover[vertex].size(); + for (unsigned int i = 0; i < sz; i++) { + simplx.push_back(cover[vertex][i]); + } + } + + std::sort(simplx.begin(), simplx.end()); + std::vector::iterator it = std::unique(simplx.begin(), simplx.end()); + simplx.resize(std::distance(simplx.begin(), it)); + simplices.push_back(simplx); + } + } + std::vector >::iterator it; + std::sort(simplices.begin(), simplices.end()); + it = std::unique(simplices.begin(), simplices.end()); + simplices.resize(std::distance(simplices.begin(), it)); + } + } + } }; -} // namespace cover_complex +} // namespace cover_complex -} // namespace Gudhi +} // namespace Gudhi #endif // GIC_H_ -- cgit v1.2.3