From 3a2f852e87514edff2d6c91c0ced2f6c53b541a1 Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Wed, 18 Mar 2020 23:43:17 +0100 Subject: [skip ci] some cleanups --- src/Collapse/example/rips_persistence_with_sc.cpp | 88 ++-------- src/Collapse/include/gudhi/PointSetGen.h | 193 ---------------------- 2 files changed, 13 insertions(+), 268 deletions(-) (limited to 'src/Collapse') diff --git a/src/Collapse/example/rips_persistence_with_sc.cpp b/src/Collapse/example/rips_persistence_with_sc.cpp index ed978daa..70f7abbb 100644 --- a/src/Collapse/example/rips_persistence_with_sc.cpp +++ b/src/Collapse/example/rips_persistence_with_sc.cpp @@ -18,34 +18,7 @@ using Rips_edge_list = Gudhi::rips_edge_list::Rips_edge_list; using Distance_matrix = std::vector>; -class extract_sub_one_skeleton -{ -public: - template - extract_sub_one_skeleton(double threshold, Filtered_sorted_edge_list & current_edge_t, Filtered_sorted_edge_list & edge_t, Fil_vector & edge_filt ) { - - auto end_it = std::upper_bound(edge_filt.begin(), edge_filt.end(), threshold); // find_index(edge_t, threshold, 0, end_idx); - size_t end_idx = std::distance(edge_filt.begin(), end_it); - for( size_t idx = 0; idx < end_idx ; idx++) { - current_edge_t.push_back(*edge_t.begin()); - edge_filt.erase(edge_filt.begin()); - edge_t.erase(edge_t.begin()); - } - - } -}; -class extract_one_new_edge -{ -public: - template - extract_one_new_edge(Filtered_sorted_edge_list & current_edge_t, Filtered_sorted_edge_list & edge_t, Fil_vector & edge_filt ) { - current_edge_t.push_back(*edge_t.begin()); - edge_filt.erase(edge_filt.begin()); - edge_t.erase(edge_t.begin()); - - } -}; class filt_edge_to_dist_matrix { @@ -94,53 +67,42 @@ int main(int argc, char * const argv[]) { char manifold; Vector_of_points * point_vector; - Vector_of_points file_all_points; std::string manifold_full = "sphere"; double r_min = 0.6; int dim_max = 2; - point_generator.program_options(argc, argv, number_of_points, steps, end_threshold, repetetions, manifold, dimension, dim_max, in_file_name, out_file_name); + point_generator.program_options(argc, argv, steps, end_threshold, repetetions, manifold, dimension, dim_max, in_file_name, out_file_name); std::cout << "The current input values to run the program is: "<< std::endl; - std::cout << "number_of_points, steps, end_threshold, repetetions, manifold, dimension, max_complex_dimension, in_file_name, out_file_name" << std::endl; - std::cout << number_of_points << ", " << steps << ", " << end_threshold << ", " << repetetions << ", " << manifold << ", " << dimension << ", " << dim_max << ", " << in_file_name << ", " << out_file_name << std::endl; + std::cout << "steps, end_threshold, repetetions, manifold, dimension, max_complex_dimension, in_file_name, out_file_name" << std::endl; + std::cout << steps << ", " << end_threshold << ", " << repetetions << ", " << manifold << ", " << dimension << ", " << dim_max << ", " << in_file_name << ", " << out_file_name << std::endl; - if(manifold == 'f' || manifold =='F') { - Gudhi::Points_off_reader off_reader(in_file_name); - if (!off_reader.is_valid()) { - std::cerr << "Unable to read file " << in_file_name << "\n"; - exit(-1); // ----- >> - } - - file_all_points = Vector_of_points(off_reader.get_point_cloud()); - dimension = file_all_points[0].dimension() ; - std::cout << "Successfully read " << file_all_points.size() << " point_vector.\n"; - std::cout << "Ambient dimension is " << dimension << ".\n"; - } - Map map_empty; std::string filediag_aft ("./PersistenceOutput/collapsed_persistence_diags") ; - std::string origPoints ("./PersistenceOutput/pointsamaple.off"); - // std::string otherStats ("./PersistenceOutput/maximal_simplx_cnt"); - // otherStats = otherStats+"_"+ out_file_name+ ".txt"; filediag_aft = filediag_aft+"_"+ out_file_name+ ".txt"; double currentCreationTime = 0.0; - point_vector = new Vector_of_points(); Distance_matrix distances; Distance_matrix *sparse_distances = new Distance_matrix(); if(manifold == 'f' || manifold =='f') { - // Subsampling from all points for each iterations - Gudhi::subsampling::pick_n_random_points(file_all_points, number_of_points, std::back_inserter(*point_vector)); + Gudhi::Points_off_reader off_reader(in_file_name); + if (!off_reader.is_valid()) { + std::cerr << "Unable to read file " << in_file_name << "\n"; + exit(-1); // ----- >> + } + + point_vector = new Vector_of_points(off_reader.get_point_cloud().begin(), off_reader.get_point_cloud().end()); + dimension = point_vector->at(0).dimension() ; number_of_points = point_vector->size(); - std::cout << number_of_points << " points succesfully chosen randomly of dimension "<< dimension << " ." << std::endl; + std::cout << "Successfully read " << number_of_points << " point_vector.\n"; + std::cout << "Ambient dimension is " << dimension << ".\n"; } else if (manifold == 'm'){ std::string csv_file_name(in_file_name); @@ -155,10 +117,6 @@ int main(int argc, char * const argv[]) { std::cout << "Point Set Generated." <at(i)); - // point_generator.output_points(*point_vector, origPoints); - Filtered_sorted_edge_list * edge_t = new Filtered_sorted_edge_list(); std::cout << "Computing the one-skeleton for threshold: " << end_threshold << std::endl; @@ -167,22 +125,12 @@ int main(int argc, char * const argv[]) { Rips_edge_list Rips_edge_list_from_file(distances, end_threshold); Rips_edge_list_from_file.create_edges(*edge_t); std::cout<< "Sorted edge list computed" << std::endl; - - //Creating the Rips Complex - //Rips_complex rips_complex_from_file(distances, end_threshold); - //rips_complex_from_file.create_complex(*subComplex, dim_max); - //std::cout<< "Rips complex computed" << std::endl; } else{ //Point cloud input //Creating the edge list Rips_edge_list Rips_edge_list_from_points(*point_vector, end_threshold, Gudhi::Euclidean_distance()); Rips_edge_list_from_points.create_edges(*edge_t); std::cout<< "Sorted edge list computed" << std::endl; std::cout << "Total number of edges before collapse are: " << edge_t->size() << std::endl; - - // Creating the Rips Complex - // Rips_complex rips_complex_after_collapse(*point_vector, end_threshold, Gudhi::Euclidean_distance()); - // rips_complex_from_points.create_complex(*subComplex, dim_max); - // std::cout<< "Rips complex computed" << std::endl; } //Now we will perform filtered edge collapse to sparsify the edge list edge_t. @@ -207,28 +155,18 @@ int main(int argc, char * const argv[]) { // Rips_complex rips_complex_before_collapse(distances, end_threshold); Rips_complex rips_complex_after_collapse(*sparse_distances, end_threshold); - // Rips_complex rips_complex_after_collapse(*point_vector, end_threshold, Gudhi::Euclidean_distance()); - // Construct the Rips complex in a Simplex Tree Simplex_tree simplex_tree_aft; - // Simplex_tree simplex_tree_bfr; - // rips_complex_before_collapse.create_complex(simplex_tree_bfr, dim_max); rips_complex_after_collapse.create_complex(simplex_tree_aft, dim_max); - // std::cout << "The complex contains " << simplex_tree_bfr.num_simplices() << " simplices before collapse. \n"; - // std::cout << " and has dimension " << simplex_tree_bfr.dimension() << " \n"; - std::cout << "The complex contains " << simplex_tree_aft.num_simplices() << " simplices after collapse. \n"; std::cout << " and has dimension " << simplex_tree_aft.dimension() << " \n"; // Sort the simplices in the order of the filtration - // simplex_tree_bfr.initialize_filtration(); simplex_tree_aft.initialize_filtration(); // Compute the persistence diagram of the complex - // Persistent_cohomology pcoh_bfr(simplex_tree_bfr); Persistent_cohomology pcoh_aft(simplex_tree_aft); // initializes the coefficient field for homology - // pcoh_bfr.init_coefficients(2); pcoh_aft.init_coefficients(2); // pcoh_bfr.compute_persistent_cohomology(steps); diff --git a/src/Collapse/include/gudhi/PointSetGen.h b/src/Collapse/include/gudhi/PointSetGen.h index a7d6956d..be83090b 100644 --- a/src/Collapse/include/gudhi/PointSetGen.h +++ b/src/Collapse/include/gudhi/PointSetGen.h @@ -14,11 +14,8 @@ const double PI = 3.141592653589793238463; #define _USE_MATH_DEFINES class PointSetGen { - private: - double unirand(){return (double) rand()/(double) RAND_MAX;} public: void program_options(int argc, char * const argv[] - , std::size_t & number_of_points , double & steps , double & end_thresold , int & repetetions @@ -33,8 +30,6 @@ class PointSetGen { po::options_description visible("Allowed options", 100); visible.add_options() ("help,h", "produce help message") - ("number,n", po::value(&number_of_points)->default_value(0), - "Number of generated point_vector.") ("steps,s", po::value(&steps)->default_value(0.1), "Steps of the threshold") ("end_thresold,e", po::value(&end_thresold)->default_value(1), @@ -78,192 +73,4 @@ class PointSetGen { std::abort(); } } - - - void generate_points_sphere(Vector_of_points& W, int nbP, int dim, double radius) { - CGAL::Random_points_on_sphere_d rp(dim+1, radius); - for (int i = 0; i < nbP; i++) - W.push_back(*rp++); - } - // void generate_fibonaci_grid_sphere(Vector_of_points& W, int nbP, int dim, double radius) - // { - - // } - - void generate_grid_2sphere(Vector_of_points& W, int nbP, int r ) - { - std::vector coords; - int Ncount = 0; - double p,v; //the angles phi and psi - int M_p; - - double a = (4*PI*pow(r,2))/nbP; - double d = sqrt(a); - - int M_v = PI/d; - - double d_v = PI/M_v; - double d_p = a/d_v; - - for( int m = 0; m < M_v ; m++) { - v = (PI*(m + 0.5))/M_v; - M_p = ((2*PI*sin(v))/d_p); - for(int n = 0; n < M_p ; n++) { - p = (2*PI*n)/M_p; - coords = {r*sin(v)*cos(p), r*sin(v)*sin(p), r*cos(v)}; - W.push_back(Point(coords)); - Ncount += 1; - } - } - } - - - /* - Generates point sets on spheres wedged at origin, sphere can have different radii from with steps of - Number of points on the sphere can also be different from for the smallest sphere and then multiplied by - */ - void generate_points_wedged_sphere(Vector_of_points& W, int init_nbP, int dim, double init_radius, int multiplier_step, int nbSpheres) { - double radius = init_radius; - int nbP = init_nbP; - std::vector translation; - for(int d = 0; d< dim; d++) { - translation.push_back(0); - } - for(int s = 0; s < nbSpheres; s++) { - CGAL::Random_points_on_sphere_d rp(dim+1, radius); - for (int i = 0; i < nbP; i++) { - W.push_back(add_point(*rp++, translation, dim)); - } - nbP = nbP*multiplier_step; - radius = radius*multiplier_step; - translation.at(dim-1) = (radius - init_radius); - } - } - - void generate_points_concentric_sphere(Vector_of_points& W, int init_nbP, int dim, int init_radius, int multiplier_step, int nbSpheres) { - double radius = init_radius; - int nbP = init_nbP; - - for(int s = 0; s < nbSpheres; s++) { - CGAL::Random_points_on_sphere_d rp(dim+1, radius); - for (int i = 0; i < nbP; i++) { - W.push_back(*rp++); - } - nbP = nbP*(pow(multiplier_step,2)); - radius = radius*multiplier_step; - } - - } - void generate_points_2annulus(Vector_of_points& W, int nbP, double r_min, double r_max) { - double rho, theta; - double x, y; - std::vector coords; - double r_min_sq = pow(r_min,2); - double r_max_sq = pow(r_max,2); - - srand(time(NULL)); - for (int i=0; i coords; - double r_min_cube = pow(r_min,3); - double r_max_cube = pow(r_max,3); - - srand(time(NULL)); - for (int i=0; i rp(dim, radius); - for (int i = 0; i < nbP; i++) - W.push_back(*rp++); - } - - void generate_points_cube(Vector_of_points& W, int nbP, int dim) { - CGAL::Random_points_in_cube_d rp(dim, 6); - for (int i = 0; i < nbP; i++) - W.push_back(*rp++); - } - - void add_point_vectors(Vector_of_points& V, Vector_of_points& U, int nbP, int dim) { // Adds two point vectors of the same size (nbP), by Modifying the first one, V = V+W. - for (int i = 0; i < nbP; i++) - { - V[i] = add_point(V[i], U[i], dim); - } - } - - //returns x = x+y; - Point add_point(const Point & x, const Point & y, int dim) { - std::vector coords; - for(int i =0; i< dim; i++) - coords.push_back(x[i]+y[i]); - return Point(coords); - } - void print_point(const Point & x) { - std::cout<< "("; - for(auto & p : x){ - std::cout<< p << ", " ; - } - std::cout<< ")" << std::endl; - } - void output_points(Vector_of_points & W, std::string outFile) { - std::ofstream myfile (outFile, std::ios::app); - if (myfile.is_open()) { - myfile << "OFF" << " " << W.size() << " " << W.at(0).size() << std::endl; - for(auto & v : W){ - for(auto & x : v){ - myfile<< x << " " ; - } - myfile << std::endl; - } - myfile << "# Tower updated for the additional subcomplex.\n"; - myfile.close(); - } - else { - std::cerr << "Unable to open file"; - exit(-1) ; - } - } - - Point noise_point(double noise_param, int dim, double radius) { - std::vector noise; - for(int d = 0; d< dim; d++){ - if(d % 2) - noise.push_back(-noise_param*radius); - else - noise.push_back(noise_param*radius); - } - return Point(noise); - } - //add noise to the points in W. - void add_noise(Vector_of_points& W, int nbP, int dim, double radius, double noise_param) { - Point noise = noise_point(noise_param, dim, radius); - for(Vector_of_points::iterator it = W.begin(); it != W.end(); it++ ) { - *it = add_point(*it,noise,dim); - } - } - - PointSetGen(){} - ~PointSetGen(){} }; \ No newline at end of file -- cgit v1.2.3