diff options
Diffstat (limited to 'src/Witness_complex/example')
6 files changed, 421 insertions, 0 deletions
diff --git a/src/Witness_complex/example/CMakeLists.txt b/src/Witness_complex/example/CMakeLists.txt new file mode 100644 index 00000000..5860f3a3 --- /dev/null +++ b/src/Witness_complex/example/CMakeLists.txt @@ -0,0 +1,34 @@ +project(Witness_complex_examples) + +add_executable ( Witness_complex_example_nearest_landmark_table example_nearest_landmark_table.cpp ) +if (TBB_FOUND) + target_link_libraries(Witness_complex_example_nearest_landmark_table ${TBB_LIBRARIES}) +endif() +add_test(NAME Witness_complex_example_nearest_landmark_table + COMMAND $<TARGET_FILE:Witness_complex_example_nearest_landmark_table>) + +install(TARGETS Witness_complex_example_nearest_landmark_table DESTINATION bin) + +# CGAL and Eigen3 are required for Euclidean version of Witness +if(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + add_executable( Witness_complex_example_off example_witness_complex_off.cpp ) + add_executable ( Witness_complex_example_sphere example_witness_complex_sphere.cpp ) + + add_executable( Witness_complex_example_strong_off example_strong_witness_complex_off.cpp ) + target_link_libraries(Witness_complex_example_strong_off) + + add_test(NAME Witness_complex_example_off_test_torus + COMMAND $<TARGET_FILE:Witness_complex_example_off> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "20" "1.0" "3") + add_test(NAME Witness_complex_example_test_sphere_10 + COMMAND $<TARGET_FILE:Witness_complex_example_sphere> "10") + add_test(NAME Witness_complex_example_strong_off_test_torus + COMMAND $<TARGET_FILE:Witness_complex_example_strong_off> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "20" "1.0" "3") + + install(TARGETS Witness_complex_example_off DESTINATION bin) + install(TARGETS Witness_complex_example_sphere DESTINATION bin) + install(TARGETS Witness_complex_example_strong_off DESTINATION bin) + + +endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) diff --git a/src/Witness_complex/example/example_nearest_landmark_table.cpp b/src/Witness_complex/example/example_nearest_landmark_table.cpp new file mode 100644 index 00000000..441900c1 --- /dev/null +++ b/src/Witness_complex/example/example_nearest_landmark_table.cpp @@ -0,0 +1,44 @@ +#define BOOST_PARAMETER_MAX_ARITY 12 + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Witness_complex.h> +#include <gudhi/Persistent_cohomology.h> + +#include <iostream> +#include <fstream> +#include <utility> +#include <string> +#include <vector> + +int main(int argc, char * const argv[]) { + using Nearest_landmark_range = std::vector<std::pair<std::size_t, double>>; + using Nearest_landmark_table = std::vector<Nearest_landmark_range>; + using Witness_complex = Gudhi::witness_complex::Witness_complex<Nearest_landmark_table>; + using Simplex_tree = Gudhi::Simplex_tree<>; + using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; + using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>; + + Simplex_tree simplex_tree; + + // Example contains 5 witnesses and 5 landmarks + Nearest_landmark_table nlt = { + {{0, 0.0}, {1, 0.1}, {2, 0.2}, {3, 0.3}, {4, 0.4}}, // witness 0 + {{1, 0.0}, {2, 0.1}, {3, 0.2}, {4, 0.3}, {0, 0.4}}, // witness 1 + {{2, 0.0}, {3, 0.1}, {4, 0.2}, {0, 0.3}, {1, 0.4}}, // witness 2 + {{3, 0.0}, {4, 0.1}, {0, 0.2}, {1, 0.3}, {2, 0.4}}, // witness 3 + {{4, 0.0}, {0, 0.1}, {1, 0.2}, {2, 0.3}, {3, 0.4}} // witness 4 + }; + /* distance(witness3, landmark3) is 0, distance(witness3, landmark4) is 0.1, etc. */ + + Witness_complex witness_complex(nlt); + witness_complex.create_complex(simplex_tree, .41); + + std::cout << "Number of simplices: " << simplex_tree.num_simplices() << std::endl; + + Persistent_cohomology pcoh(simplex_tree); + // initializes the coefficient field for homology + pcoh.init_coefficients(11); + + pcoh.compute_persistent_cohomology(-0.1); + pcoh.output_diagram(); +} diff --git a/src/Witness_complex/example/example_strong_witness_complex_off.cpp b/src/Witness_complex/example/example_strong_witness_complex_off.cpp new file mode 100644 index 00000000..19f73836 --- /dev/null +++ b/src/Witness_complex/example/example_strong_witness_complex_off.cpp @@ -0,0 +1,57 @@ +#include <gudhi/Simplex_tree.h> +#include <gudhi/Euclidean_strong_witness_complex.h> +#include <gudhi/pick_n_random_points.h> +#include <gudhi/choose_n_farthest_points.h> +#include <gudhi/Points_off_io.h> + +#include <CGAL/Epick_d.h> + +#include <iostream> +#include <fstream> +#include <ctime> +#include <string> +#include <vector> + +using K = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>; +using Point_d = typename K::Point_d; +using Witness_complex = Gudhi::witness_complex::Euclidean_strong_witness_complex<K>; +using Point_vector = std::vector<Point_d>; + +int main(int argc, char* const argv[]) { + if (argc != 5) { + std::cerr << "Usage: " << argv[0] << " path_to_point_file number_of_landmarks max_squared_alpha limit_dimension\n"; + return 0; + } + + std::string file_name = argv[1]; + int nbL = atoi(argv[2]), lim_dim = atoi(argv[4]); + double alpha2 = atof(argv[3]); + clock_t start, end; + Gudhi::Simplex_tree<> simplex_tree; + + // Read the point file + Point_vector point_vector, landmarks; + Gudhi::Points_off_reader<Point_d> off_reader(file_name); + if (!off_reader.is_valid()) { + std::cerr << "Strong witness complex - Unable to read file " << file_name << "\n"; + exit(-1); // ----- >> + } + point_vector = Point_vector(off_reader.get_point_cloud()); + + std::cout << "Successfully read " << point_vector.size() << " points.\n"; + std::cout << "Ambient dimension is " << point_vector[0].dimension() << ".\n"; + + // Choose landmarks (decomment one of the following two lines) + // Gudhi::subsampling::pick_n_random_points(point_vector, nbL, std::back_inserter(landmarks)); + Gudhi::subsampling::choose_n_farthest_points(K(), point_vector, nbL, Gudhi::subsampling::random_starting_point, + std::back_inserter(landmarks)); + + // Compute witness complex + start = clock(); + Witness_complex witness_complex(landmarks, point_vector); + + witness_complex.create_complex(simplex_tree, alpha2, lim_dim); + end = clock(); + std::cout << "Strong witness complex took " << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n"; + std::cout << "Number of simplices is: " << simplex_tree.num_simplices() << "\n"; +} diff --git a/src/Witness_complex/example/example_witness_complex_off.cpp b/src/Witness_complex/example/example_witness_complex_off.cpp new file mode 100644 index 00000000..be11c955 --- /dev/null +++ b/src/Witness_complex/example/example_witness_complex_off.cpp @@ -0,0 +1,62 @@ +#include <sys/types.h> +#include <sys/stat.h> + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Euclidean_witness_complex.h> +#include <gudhi/pick_n_random_points.h> +#include <gudhi/choose_n_farthest_points.h> +#include <gudhi/Points_off_io.h> + +#include <CGAL/Epick_d.h> + +#include <iostream> +#include <fstream> +#include <ctime> +#include <string> +#include <vector> + +using K = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>; +using Point_d = K::Point_d; +using Witness_complex = Gudhi::witness_complex::Euclidean_witness_complex<K>; +using Point_vector = std::vector< Point_d >; + +int main(int argc, char * const argv[]) { + if (argc != 5) { + std::cerr << "Usage: " << argv[0] + << " path_to_point_file number_of_landmarks max_squared_alpha limit_dimension\n"; + return 0; + } + + std::string file_name = argv[1]; + int nbL = atoi(argv[2]), lim_dim = atoi(argv[4]); + double alpha2 = atof(argv[3]); + clock_t start, end; + Gudhi::Simplex_tree<> simplex_tree; + + // Read the point file + Point_vector point_vector, landmarks; + Gudhi::Points_off_reader<Point_d> off_reader(file_name); + if (!off_reader.is_valid()) { + std::cerr << "Witness complex - Unable to read file " << file_name << "\n"; + exit(-1); // ----- >> + } + point_vector = Point_vector(off_reader.get_point_cloud()); + + std::cout << "Successfully read " << point_vector.size() << " points.\n"; + std::cout << "Ambient dimension is " << point_vector[0].dimension() << ".\n"; + + // Choose landmarks (decomment one of the following two lines) + // Gudhi::subsampling::pick_n_random_points(point_vector, nbL, std::back_inserter(landmarks)); + Gudhi::subsampling::choose_n_farthest_points(K(), point_vector, nbL, Gudhi::subsampling::random_starting_point, std::back_inserter(landmarks)); + + // Compute witness complex + start = clock(); + Witness_complex witness_complex(landmarks, + point_vector); + + witness_complex.create_complex(simplex_tree, alpha2, lim_dim); + end = clock(); + std::cout << "Witness complex took " + << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n"; + std::cout << "Number of simplices is: " << simplex_tree.num_simplices() << "\n"; +} diff --git a/src/Witness_complex/example/example_witness_complex_sphere.cpp b/src/Witness_complex/example/example_witness_complex_sphere.cpp new file mode 100644 index 00000000..9e3c972d --- /dev/null +++ b/src/Witness_complex/example/example_witness_complex_sphere.cpp @@ -0,0 +1,70 @@ +#define BOOST_PARAMETER_MAX_ARITY 12 + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Euclidean_witness_complex.h> +#include <gudhi/pick_n_random_points.h> +#include <gudhi/choose_n_farthest_points.h> +#include <gudhi/reader_utils.h> + +#include <CGAL/Epick_d.h> + +#include <iostream> +#include <fstream> +#include <ctime> +#include <utility> +#include <string> +#include <vector> + +#include "generators.h" + +/** Write a gnuplot readable file. + * Data range is a random access range of pairs (arg, value) + */ +template <typename Data_range> +void write_data(Data_range& data, std::string filename) { + std::ofstream ofs(filename, std::ofstream::out); + for (auto entry : data) ofs << entry.first << ", " << entry.second << "\n"; + ofs.close(); +} + +int main(int argc, char* const argv[]) { + using Kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>; + using Witness_complex = Gudhi::witness_complex::Euclidean_witness_complex<Kernel>; + + if (argc != 2) { + std::cerr << "Usage: " << argv[0] << " number_of_landmarks \n"; + return 0; + } + + int number_of_landmarks = atoi(argv[1]); + + std::vector<std::pair<int, double> > l_time; + + // Generate points + for (int nbP = 500; nbP < 10000; nbP += 500) { + clock_t start, end; + // Construct the Simplex Tree + Gudhi::Simplex_tree<> simplex_tree; + Point_Vector point_vector, landmarks; + generate_points_sphere(point_vector, nbP, 4); + std::cout << "Successfully generated " << point_vector.size() << " points.\n"; + std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; + + // Choose landmarks + start = clock(); + // Gudhi::subsampling::pick_n_random_points(point_vector, number_of_landmarks, std::back_inserter(landmarks)); + Gudhi::subsampling::choose_n_farthest_points(K(), point_vector, number_of_landmarks, + Gudhi::subsampling::random_starting_point, + std::back_inserter(landmarks)); + + // Compute witness complex + Witness_complex witness_complex(landmarks, point_vector); + witness_complex.create_complex(simplex_tree, 0); + end = clock(); + double time = static_cast<double>(end - start) / CLOCKS_PER_SEC; + std::cout << "Witness complex for " << number_of_landmarks << " landmarks took " << time << " s. \n"; + std::cout << "Number of simplices is: " << simplex_tree.num_simplices() << "\n"; + l_time.push_back(std::make_pair(nbP, time)); + } + write_data(l_time, "w_time.dat"); +} diff --git a/src/Witness_complex/example/generators.h b/src/Witness_complex/example/generators.h new file mode 100644 index 00000000..1900e1e4 --- /dev/null +++ b/src/Witness_complex/example/generators.h @@ -0,0 +1,154 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Siargey Kachanovich + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef GENERATORS_H_ +#define GENERATORS_H_ + +#include <CGAL/Epick_d.h> +#include <CGAL/point_generators_d.h> +#include <CGAL/Random.h> + +#include <fstream> +#include <string> +#include <vector> +#include <cmath> + +using K = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>; +using FT = K::FT; +using Point_d = K::Point_d; +using Point_Vector = std::vector<Point_d>; +using Random_cube_iterator = CGAL::Random_points_in_cube_d<Point_d>; +using Random_point_iterator = CGAL::Random_points_in_ball_d<Point_d>; + +/** + * \brief Rock age method of reading off file + * + */ +inline void +off_reader_cust(std::string file_name, std::vector<Point_d> & points) { + std::ifstream in_file(file_name.c_str(), std::ios::in); + if (!in_file.is_open()) { + std::cerr << "Unable to open file " << file_name << std::endl; + return; + } + std::string line; + double x; + // Line OFF. No need in it + if (!getline(in_file, line)) { + std::cerr << "No line OFF\n"; + return; + } + // Line with 3 numbers. No need + if (!getline(in_file, line)) { + std::cerr << "No line with 3 numbers\n"; + return; + } + // Reading points + while (getline(in_file, line)) { + std::vector< double > point; + std::istringstream iss(line); + while (iss >> x) { + point.push_back(x); + } + points.push_back(Point_d(point)); + } + in_file.close(); +} + +/** + * \brief Customized version of read_points + * which takes into account a possible nbP first line + * + */ +inline void +read_points_cust(std::string file_name, Point_Vector & points) { + std::ifstream in_file(file_name.c_str(), std::ios::in); + if (!in_file.is_open()) { + std::cerr << "Unable to open file " << file_name << std::endl; + return; + } + std::string line; + double x; + while (getline(in_file, line)) { + std::vector< double > point; + std::istringstream iss(line); + while (iss >> x) { + point.push_back(x); + } + Point_d p(point.begin(), point.end()); + if (point.size() != 1) + points.push_back(p); + } + in_file.close(); +} + +/** \brief Generate points on a grid in a cube of side 2 + * having {+-1}^D as vertices and insert them in W. + * The grid has "width" points on each side. + * If torus is true then it is supposed that the cube represents + * a flat torus, hence the opposite borders are associated. + * The points on border in this case are not placed twice. + */ +void generate_points_grid(Point_Vector& W, int width, int D, bool torus) { + int nb_points = 1; + for (int i = 0; i < D; ++i) + nb_points *= width; + for (int i = 0; i < nb_points; ++i) { + std::vector<double> point; + int cell_i = i; + for (int l = 0; l < D; ++l) { + if (torus) + point.push_back(-1 + (2.0 / (width - 1))*(cell_i % width)); + else + point.push_back(-1 + (2.0 / width)*(cell_i % width)); + // attention: the bottom and the right are covered too! + cell_i /= width; + } + W.push_back(point); + } +} + +/** \brief Generate nbP points uniformly in a cube of side 2 + * having {+-1}^dim as its vertices and insert them in W. + */ +void generate_points_random_box(Point_Vector& W, int nbP, int dim) { + Random_cube_iterator rp(dim, 1.0); + for (int i = 0; i < nbP; i++) { + W.push_back(*rp++); + } +} + +/** \brief Generate nbP points uniformly on a (dim-1)-sphere + * and insert them in W. + */ +void generate_points_sphere(Point_Vector& W, int nbP, int dim) { + CGAL::Random_points_on_sphere_d<Point_d> rp(dim, 1); + for (int i = 0; i < nbP; i++) + W.push_back(*rp++); +} + +/** \brief Generate nbP points on a (flat) d-torus embedded in R^{2d} + * + */ +void generate_points_torus(Point_Vector& W, int nbP, int dim) { + CGAL::Random rand; + const double pi = std::acos(-1); + for (int i = 0; i < nbP; i++) { + std::vector<FT> point; + for (int j = 0; j < dim; j++) { + double alpha = rand.uniform_real(static_cast<double>(0), 2*pi); + point.push_back(sin(alpha)); + point.push_back(cos(alpha)); + } + W.push_back(Point_d(point)); + } +} + +#endif // GENERATORS_H_ |