diff options
Diffstat (limited to 'src/common/include/gudhi')
-rw-r--r-- | src/common/include/gudhi/Clock.h | 48 | ||||
-rw-r--r-- | src/common/include/gudhi/Debug_utils.h | 2 | ||||
-rw-r--r-- | src/common/include/gudhi/console_color.h | 97 | ||||
-rw-r--r-- | src/common/include/gudhi/distance_functions.h | 36 | ||||
-rw-r--r-- | src/common/include/gudhi/graph_simplicial_complex.h | 59 | ||||
-rw-r--r-- | src/common/include/gudhi/random_point_generators.h | 474 | ||||
-rw-r--r-- | src/common/include/gudhi/reader_utils.h | 166 |
7 files changed, 758 insertions, 124 deletions
diff --git a/src/common/include/gudhi/Clock.h b/src/common/include/gudhi/Clock.h index 04c6ffb9..77f196ca 100644 --- a/src/common/include/gudhi/Clock.h +++ b/src/common/include/gudhi/Clock.h @@ -27,47 +27,55 @@ #include <string> +namespace Gudhi { + class Clock { public: - Clock() : end_called(false) { - startTime = boost::posix_time::microsec_clock::local_time(); - } - - Clock(const std::string& msg_) { - end_called = false; - begin(); - msg = msg_; - } + // Construct and start the timer + Clock(const std::string& msg_ = std::string()) + : startTime(boost::posix_time::microsec_clock::local_time()), + end_called(false), + msg(msg_) { } + // Restart the timer void begin() const { end_called = false; startTime = boost::posix_time::microsec_clock::local_time(); } + // Stop the timer void end() const { end_called = true; endTime = boost::posix_time::microsec_clock::local_time(); } + std::string message() const { + return msg; + } + + // Print current value to std::cout void print() const { std::cout << *this << std::endl; } friend std::ostream& operator<<(std::ostream& stream, const Clock& clock) { - if (!clock.end_called) - clock.end(); + if (!clock.msg.empty()) + stream << clock.msg << ": "; - if (!clock.end_called) { - stream << "end not called"; - } else { - stream << clock.msg << ":" << clock.num_seconds() << "s"; - } + stream << clock.num_seconds() << "s"; return stream; } + // Get the number of seconds between the timer start and: + // - the last call of end() if it was called + // - or now otherwise. In this case, the timer is not stopped. double num_seconds() const { - if (!end_called) return -1; - return (endTime - startTime).total_milliseconds() / 1000.; + if (!end_called) { + auto end = boost::posix_time::microsec_clock::local_time(); + return (end - startTime).total_milliseconds() / 1000.; + } else { + return (endTime - startTime).total_milliseconds() / 1000.; + } } private: @@ -76,4 +84,6 @@ class Clock { std::string msg; }; -#endif // CLOCK_H_ +} // namespace Gudhi + +#endif // CLOCK_H_ diff --git a/src/common/include/gudhi/Debug_utils.h b/src/common/include/gudhi/Debug_utils.h index 7573a9db..8ed3b7b3 100644 --- a/src/common/include/gudhi/Debug_utils.h +++ b/src/common/include/gudhi/Debug_utils.h @@ -33,8 +33,10 @@ // Could assert in release mode, but cmake sets NDEBUG (for "NO DEBUG") in this mode, means assert does nothing. #ifdef GUDHI_DEBUG #define GUDHI_CHECK(expression, excpt) if ((expression) == 0) throw excpt + #define GUDHI_CHECK_code(CODE) CODE #else #define GUDHI_CHECK(expression, excpt) (void) 0 + #define GUDHI_CHECK_code(CODE) #endif #define PRINT(a) std::cerr << #a << ": " << (a) << " (DISP)" << std::endl diff --git a/src/common/include/gudhi/console_color.h b/src/common/include/gudhi/console_color.h new file mode 100644 index 00000000..c4671da3 --- /dev/null +++ b/src/common/include/gudhi/console_color.h @@ -0,0 +1,97 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Clement Jamin + * + * Copyright (C) 2016 INRIA Sophia-Antipolis (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef CONSOLE_COLOR_H_ +#define CONSOLE_COLOR_H_ + +#include <iostream> + +#if defined(WIN32) +#include <windows.h> +#endif + +inline std::ostream& blue(std::ostream &s) { +#if defined(WIN32) + HANDLE hStdout = GetStdHandle(STD_OUTPUT_HANDLE); + SetConsoleTextAttribute(hStdout, + FOREGROUND_BLUE | FOREGROUND_GREEN | FOREGROUND_INTENSITY); +#else + s << "\x1b[0;34m"; +#endif + return s; +} + +inline std::ostream& red(std::ostream &s) { +#if defined(WIN32) + HANDLE hStdout = GetStdHandle(STD_OUTPUT_HANDLE); + SetConsoleTextAttribute(hStdout, FOREGROUND_RED | FOREGROUND_INTENSITY); +#else + s << "\x1b[0;31m"; +#endif + return s; +} + +inline std::ostream& green(std::ostream &s) { +#if defined(WIN32) + HANDLE hStdout = GetStdHandle(STD_OUTPUT_HANDLE); + SetConsoleTextAttribute(hStdout, FOREGROUND_GREEN | FOREGROUND_INTENSITY); +#else + s << "\x1b[0;32m"; +#endif + return s; +} + +inline std::ostream& yellow(std::ostream &s) { +#if defined(WIN32) + HANDLE hStdout = GetStdHandle(STD_OUTPUT_HANDLE); + SetConsoleTextAttribute(hStdout, + FOREGROUND_GREEN | FOREGROUND_RED | FOREGROUND_INTENSITY); +#else + s << "\x1b[0;33m"; +#endif + return s; +} + +inline std::ostream& white(std::ostream &s) { +#if defined(WIN32) + HANDLE hStdout = GetStdHandle(STD_OUTPUT_HANDLE); + SetConsoleTextAttribute(hStdout, + FOREGROUND_RED | FOREGROUND_GREEN | FOREGROUND_BLUE); +#else + s << "\x1b[0;37m"; +#endif + return s; +} + +inline std::ostream& black_on_white(std::ostream &s) { +#if defined(WIN32) + HANDLE hStdout = GetStdHandle(STD_OUTPUT_HANDLE); + SetConsoleTextAttribute(hStdout, + BACKGROUND_RED | BACKGROUND_GREEN | BACKGROUND_BLUE); +#else + s << "\x1b[0;33m"; +#endif + return s; +} + + +#endif // CONSOLE_COLOR_H_ diff --git a/src/common/include/gudhi/distance_functions.h b/src/common/include/gudhi/distance_functions.h index cd518581..22747637 100644 --- a/src/common/include/gudhi/distance_functions.h +++ b/src/common/include/gudhi/distance_functions.h @@ -4,7 +4,7 @@ * * Author(s): Clément Maria * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * Copyright (C) 2014 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -24,20 +24,28 @@ #define DISTANCE_FUNCTIONS_H_ #include <cmath> // for std::sqrt +#include <type_traits> // for std::decay +#include <iterator> // for std::begin, std::end -/* Compute the Euclidean distance between two Points given - * by a range of coordinates. The points are assumed to have - * the same dimension. */ -template< typename Point > -double euclidean_distance(Point &p1, Point &p2) { - double dist = 0.; - auto it1 = p1.begin(); - auto it2 = p2.begin(); - for (; it1 != p1.end(); ++it1, ++it2) { - double tmp = *it1 - *it2; - dist += tmp*tmp; +/** @file + * @brief Global distance functions + */ + +/** @brief Compute the Euclidean distance between two Points given by a range of coordinates. The points are assumed to + * have the same dimension. */ +class Euclidean_distance { + public: + template< typename Point > + auto operator()(const Point& p1, const Point& p2) const -> typename std::decay<decltype(*std::begin(p1))>::type { + auto it1 = p1.begin(); + auto it2 = p2.begin(); + typename Point::value_type dist = 0.; + for (; it1 != p1.end(); ++it1, ++it2) { + typename Point::value_type tmp = (*it1) - (*it2); + dist += tmp*tmp; + } + return std::sqrt(dist); } - return std::sqrt(dist); -} +}; #endif // DISTANCE_FUNCTIONS_H_ diff --git a/src/common/include/gudhi/graph_simplicial_complex.h b/src/common/include/gudhi/graph_simplicial_complex.h index 042ef516..5fe7c826 100644 --- a/src/common/include/gudhi/graph_simplicial_complex.h +++ b/src/common/include/gudhi/graph_simplicial_complex.h @@ -4,7 +4,7 @@ * * Author(s): Clément Maria * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * Copyright (C) 2014 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -39,61 +39,4 @@ struct vertex_filtration_t { typedef boost::vertex_property_tag kind; }; -typedef int Vertex_handle; -typedef double Filtration_value; -typedef boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS -, boost::property < vertex_filtration_t, Filtration_value > -, boost::property < edge_filtration_t, Filtration_value > -> Graph_t; -typedef std::pair< Vertex_handle, Vertex_handle > Edge_t; - -/** \brief Output the proximity graph of the points. - * - * If points contains n elements, the proximity graph is the graph - * with n vertices, and an edge [u,v] iff the distance function between - * points u and v is smaller than threshold. - * - * The type PointCloud furnishes .begin() and .end() methods, that return - * iterators with value_type Point. - */ -template< typename PointCloud -, typename Point > -Graph_t compute_proximity_graph(PointCloud &points - , Filtration_value threshold - , Filtration_value distance(Point p1, Point p2)) { - std::vector< Edge_t > edges; - std::vector< Filtration_value > edges_fil; - std::map< Vertex_handle, Filtration_value > vertices; - - Vertex_handle idx_u, idx_v; - Filtration_value fil; - idx_u = 0; - for (auto it_u = points.begin(); it_u != points.end(); ++it_u) { - idx_v = idx_u + 1; - for (auto it_v = it_u + 1; it_v != points.end(); ++it_v, ++idx_v) { - fil = distance(*it_u, *it_v); - if (fil <= threshold) { - edges.emplace_back(idx_u, idx_v); - edges_fil.push_back(fil); - } - } - ++idx_u; - } - - Graph_t skel_graph(edges.begin() - , edges.end() - , edges_fil.begin() - , idx_u); // number of points labeled from 0 to idx_u-1 - - auto vertex_prop = boost::get(vertex_filtration_t(), skel_graph); - - boost::graph_traits<Graph_t>::vertex_iterator vi, vi_end; - for (std::tie(vi, vi_end) = boost::vertices(skel_graph); - vi != vi_end; ++vi) { - boost::put(vertex_prop, *vi, 0.); - } - - return skel_graph; -} - #endif // GRAPH_SIMPLICIAL_COMPLEX_H_ diff --git a/src/common/include/gudhi/random_point_generators.h b/src/common/include/gudhi/random_point_generators.h new file mode 100644 index 00000000..2ec465ef --- /dev/null +++ b/src/common/include/gudhi/random_point_generators.h @@ -0,0 +1,474 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Clement Jamin + * + * Copyright (C) 2016 INRIA + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef RANDOM_POINT_GENERATORS_H_ +#define RANDOM_POINT_GENERATORS_H_ + +#include <CGAL/number_utils.h> +#include <CGAL/Random.h> +#include <CGAL/point_generators_d.h> + +#include <vector> // for vector<> + +namespace Gudhi { + +/////////////////////////////////////////////////////////////////////////////// +// Note: All these functions have been tested with the CGAL::Epick_d kernel +/////////////////////////////////////////////////////////////////////////////// + +// construct_point: dim 2 + +template <typename Kernel> +typename Kernel::Point_d construct_point(const Kernel &k, + typename Kernel::FT x1, typename Kernel::FT x2) { + typename Kernel::FT tab[2]; + tab[0] = x1; + tab[1] = x2; + return k.construct_point_d_object()(2, &tab[0], &tab[2]); +} + +// construct_point: dim 3 + +template <typename Kernel> +typename Kernel::Point_d construct_point(const Kernel &k, + typename Kernel::FT x1, typename Kernel::FT x2, typename Kernel::FT x3) { + typename Kernel::FT tab[3]; + tab[0] = x1; + tab[1] = x2; + tab[2] = x3; + return k.construct_point_d_object()(3, &tab[0], &tab[3]); +} + +// construct_point: dim 4 + +template <typename Kernel> +typename Kernel::Point_d construct_point(const Kernel &k, + typename Kernel::FT x1, typename Kernel::FT x2, typename Kernel::FT x3, + typename Kernel::FT x4) { + typename Kernel::FT tab[4]; + tab[0] = x1; + tab[1] = x2; + tab[2] = x3; + tab[3] = x4; + return k.construct_point_d_object()(4, &tab[0], &tab[4]); +} + +// construct_point: dim 5 + +template <typename Kernel> +typename Kernel::Point_d construct_point(const Kernel &k, + typename Kernel::FT x1, typename Kernel::FT x2, typename Kernel::FT x3, + typename Kernel::FT x4, typename Kernel::FT x5) { + typename Kernel::FT tab[5]; + tab[0] = x1; + tab[1] = x2; + tab[2] = x3; + tab[3] = x4; + tab[4] = x5; + return k.construct_point_d_object()(5, &tab[0], &tab[5]); +} + +// construct_point: dim 6 + +template <typename Kernel> +typename Kernel::Point_d construct_point(const Kernel &k, + typename Kernel::FT x1, typename Kernel::FT x2, typename Kernel::FT x3, + typename Kernel::FT x4, typename Kernel::FT x5, typename Kernel::FT x6) { + typename Kernel::FT tab[6]; + tab[0] = x1; + tab[1] = x2; + tab[2] = x3; + tab[3] = x4; + tab[4] = x5; + tab[5] = x6; + return k.construct_point_d_object()(6, &tab[0], &tab[6]); +} + +template <typename Kernel> +std::vector<typename Kernel::Point_d> generate_points_on_plane(std::size_t num_points, int intrinsic_dim, + int ambient_dim, + double coord_min = -5., double coord_max = 5.) { + typedef typename Kernel::Point_d Point; + typedef typename Kernel::FT FT; + Kernel k; + CGAL::Random rng; + std::vector<Point> points; + points.reserve(num_points); + for (std::size_t i = 0; i < num_points;) { + std::vector<FT> pt(ambient_dim, FT(0)); + for (int j = 0; j < intrinsic_dim; ++j) + pt[j] = rng.get_double(coord_min, coord_max); + + Point p = k.construct_point_d_object()(ambient_dim, pt.begin(), pt.end()); + points.push_back(p); + ++i; + } + return points; +} + +template <typename Kernel> +std::vector<typename Kernel::Point_d> generate_points_on_moment_curve(std::size_t num_points, int dim, + typename Kernel::FT min_x, + typename Kernel::FT max_x) { + typedef typename Kernel::Point_d Point; + typedef typename Kernel::FT FT; + Kernel k; + CGAL::Random rng; + std::vector<Point> points; + points.reserve(num_points); + for (std::size_t i = 0; i < num_points;) { + FT x = rng.get_double(min_x, max_x); + std::vector<FT> coords; + coords.reserve(dim); + for (int p = 1; p <= dim; ++p) + coords.push_back(std::pow(CGAL::to_double(x), p)); + Point p = k.construct_point_d_object()( + dim, coords.begin(), coords.end()); + points.push_back(p); + ++i; + } + return points; +} + + +// R = big radius, r = small radius +template <typename Kernel/*, typename TC_basis*/> +std::vector<typename Kernel::Point_d> generate_points_on_torus_3D(std::size_t num_points, double R, double r, + bool uniform = false) { + typedef typename Kernel::Point_d Point; + typedef typename Kernel::FT FT; + Kernel k; + CGAL::Random rng; + + // if uniform + std::size_t num_lines = (std::size_t)sqrt(num_points); + + std::vector<Point> points; + points.reserve(num_points); + for (std::size_t i = 0; i < num_points;) { + FT u, v; + if (uniform) { + std::size_t k1 = i / num_lines; + std::size_t k2 = i % num_lines; + u = 6.2832 * k1 / num_lines; + v = 6.2832 * k2 / num_lines; + } else { + u = rng.get_double(0, 6.2832); + v = rng.get_double(0, 6.2832); + } + Point p = construct_point(k, + (R + r * std::cos(u)) * std::cos(v), + (R + r * std::cos(u)) * std::sin(v), + r * std::sin(u)); + points.push_back(p); + ++i; + } + return points; +} + +// "Private" function used by generate_points_on_torus_d +template <typename Kernel, typename OutputIterator> +static void generate_uniform_points_on_torus_d(const Kernel &k, int dim, std::size_t num_slices, + OutputIterator out, + double radius_noise_percentage = 0., + std::vector<typename Kernel::FT> current_point = std::vector<typename Kernel::FT>()) { + CGAL::Random rng; + int point_size = static_cast<int>(current_point.size()); + if (point_size == 2 * dim) { + *out++ = k.construct_point_d_object()(point_size, current_point.begin(), current_point.end()); + } else { + for (std::size_t slice_idx = 0; slice_idx < num_slices; ++slice_idx) { + double radius_noise_ratio = 1.; + if (radius_noise_percentage > 0.) { + radius_noise_ratio = rng.get_double( + (100. - radius_noise_percentage) / 100., + (100. + radius_noise_percentage) / 100.); + } + std::vector<typename Kernel::FT> cp2 = current_point; + double alpha = 6.2832 * slice_idx / num_slices; + cp2.push_back(radius_noise_ratio * std::cos(alpha)); + cp2.push_back(radius_noise_ratio * std::sin(alpha)); + generate_uniform_points_on_torus_d( + k, dim, num_slices, out, radius_noise_percentage, cp2); + } + } +} + +template <typename Kernel> +std::vector<typename Kernel::Point_d> generate_points_on_torus_d(std::size_t num_points, int dim, bool uniform = false, + double radius_noise_percentage = 0.) { + typedef typename Kernel::Point_d Point; + typedef typename Kernel::FT FT; + Kernel k; + CGAL::Random rng; + + std::vector<Point> points; + points.reserve(num_points); + if (uniform) { + std::size_t num_slices = (std::size_t)std::pow(num_points, 1. / dim); + generate_uniform_points_on_torus_d( + k, dim, num_slices, std::back_inserter(points), radius_noise_percentage); + } else { + for (std::size_t i = 0; i < num_points;) { + double radius_noise_ratio = 1.; + if (radius_noise_percentage > 0.) { + radius_noise_ratio = rng.get_double( + (100. - radius_noise_percentage) / 100., + (100. + radius_noise_percentage) / 100.); + } + std::vector<typename Kernel::FT> pt; + pt.reserve(dim * 2); + for (int curdim = 0; curdim < dim; ++curdim) { + FT alpha = rng.get_double(0, 6.2832); + pt.push_back(radius_noise_ratio * std::cos(alpha)); + pt.push_back(radius_noise_ratio * std::sin(alpha)); + } + + Point p = k.construct_point_d_object()(pt.begin(), pt.end()); + points.push_back(p); + ++i; + } + } + return points; +} + +template <typename Kernel> +std::vector<typename Kernel::Point_d> generate_points_on_sphere_d(std::size_t num_points, int dim, double radius, + double radius_noise_percentage = 0.) { + typedef typename Kernel::Point_d Point; + Kernel k; + CGAL::Random rng; + CGAL::Random_points_on_sphere_d<Point> generator(dim, radius); + std::vector<Point> points; + points.reserve(num_points); + for (std::size_t i = 0; i < num_points;) { + Point p = *generator++; + if (radius_noise_percentage > 0.) { + double radius_noise_ratio = rng.get_double( + (100. - radius_noise_percentage) / 100., + (100. + radius_noise_percentage) / 100.); + + typename Kernel::Point_to_vector_d k_pt_to_vec = + k.point_to_vector_d_object(); + typename Kernel::Vector_to_point_d k_vec_to_pt = + k.vector_to_point_d_object(); + typename Kernel::Scaled_vector_d k_scaled_vec = + k.scaled_vector_d_object(); + p = k_vec_to_pt(k_scaled_vec(k_pt_to_vec(p), radius_noise_ratio)); + } + points.push_back(p); + ++i; + } + return points; +} + +template <typename Kernel> +std::vector<typename Kernel::Point_d> generate_points_on_two_spheres_d(std::size_t num_points, int dim, double radius, + double distance_between_centers, + double radius_noise_percentage = 0.) { + typedef typename Kernel::FT FT; + typedef typename Kernel::Point_d Point; + typedef typename Kernel::Vector_d Vector; + Kernel k; + CGAL::Random rng; + CGAL::Random_points_on_sphere_d<Point> generator(dim, radius); + std::vector<Point> points; + points.reserve(num_points); + + std::vector<FT> t(dim, FT(0)); + t[0] = distance_between_centers; + Vector c1_to_c2(t.begin(), t.end()); + + for (std::size_t i = 0; i < num_points;) { + Point p = *generator++; + if (radius_noise_percentage > 0.) { + double radius_noise_ratio = rng.get_double( + (100. - radius_noise_percentage) / 100., + (100. + radius_noise_percentage) / 100.); + + typename Kernel::Point_to_vector_d k_pt_to_vec = + k.point_to_vector_d_object(); + typename Kernel::Vector_to_point_d k_vec_to_pt = + k.vector_to_point_d_object(); + typename Kernel::Scaled_vector_d k_scaled_vec = + k.scaled_vector_d_object(); + p = k_vec_to_pt(k_scaled_vec(k_pt_to_vec(p), radius_noise_ratio)); + } + + typename Kernel::Translated_point_d k_transl = + k.translated_point_d_object(); + Point p2 = k_transl(p, c1_to_c2); + points.push_back(p); + points.push_back(p2); + i += 2; + } + return points; +} + +// Product of a 3-sphere and a circle => d = 3 / D = 5 + +template <typename Kernel> +std::vector<typename Kernel::Point_d> generate_points_on_3sphere_and_circle(std::size_t num_points, + double sphere_radius) { + typedef typename Kernel::FT FT; + typedef typename Kernel::Point_d Point; + Kernel k; + CGAL::Random rng; + CGAL::Random_points_on_sphere_d<Point> generator(3, sphere_radius); + std::vector<Point> points; + points.reserve(num_points); + + typename Kernel::Compute_coordinate_d k_coord = + k.compute_coordinate_d_object(); + for (std::size_t i = 0; i < num_points;) { + Point p_sphere = *generator++; // First 3 coords + + FT alpha = rng.get_double(0, 6.2832); + std::vector<FT> pt(5); + pt[0] = k_coord(p_sphere, 0); + pt[1] = k_coord(p_sphere, 1); + pt[2] = k_coord(p_sphere, 2); + pt[3] = std::cos(alpha); + pt[4] = std::sin(alpha); + Point p(pt.begin(), pt.end()); + points.push_back(p); + ++i; + } + return points; +} + +// a = big radius, b = small radius +template <typename Kernel> +std::vector<typename Kernel::Point_d> generate_points_on_klein_bottle_3D(std::size_t num_points, double a, double b, + bool uniform = false) { + typedef typename Kernel::Point_d Point; + typedef typename Kernel::FT FT; + Kernel k; + CGAL::Random rng; + + // if uniform + std::size_t num_lines = (std::size_t)sqrt(num_points); + + std::vector<Point> points; + points.reserve(num_points); + for (std::size_t i = 0; i < num_points;) { + FT u, v; + if (uniform) { + std::size_t k1 = i / num_lines; + std::size_t k2 = i % num_lines; + u = 6.2832 * k1 / num_lines; + v = 6.2832 * k2 / num_lines; + } else { + u = rng.get_double(0, 6.2832); + v = rng.get_double(0, 6.2832); + } + double tmp = cos(u / 2) * sin(v) - sin(u / 2) * sin(2. * v); + Point p = construct_point(k, + (a + b * tmp) * cos(u), + (a + b * tmp) * sin(u), + b * (sin(u / 2) * sin(v) + cos(u / 2) * sin(2. * v))); + points.push_back(p); + ++i; + } + return points; +} + +// a = big radius, b = small radius +template <typename Kernel> +std::vector<typename Kernel::Point_d> generate_points_on_klein_bottle_4D(std::size_t num_points, double a, double b, + double noise = 0., bool uniform = false) { + typedef typename Kernel::Point_d Point; + typedef typename Kernel::FT FT; + Kernel k; + CGAL::Random rng; + + // if uniform + std::size_t num_lines = (std::size_t)sqrt(num_points); + + std::vector<Point> points; + points.reserve(num_points); + for (std::size_t i = 0; i < num_points;) { + FT u, v; + if (uniform) { + std::size_t k1 = i / num_lines; + std::size_t k2 = i % num_lines; + u = 6.2832 * k1 / num_lines; + v = 6.2832 * k2 / num_lines; + } else { + u = rng.get_double(0, 6.2832); + v = rng.get_double(0, 6.2832); + } + Point p = construct_point(k, + (a + b * cos(v)) * cos(u) + (noise == 0. ? 0. : rng.get_double(0, noise)), + (a + b * cos(v)) * sin(u) + (noise == 0. ? 0. : rng.get_double(0, noise)), + b * sin(v) * cos(u / 2) + (noise == 0. ? 0. : rng.get_double(0, noise)), + b * sin(v) * sin(u / 2) + (noise == 0. ? 0. : rng.get_double(0, noise))); + points.push_back(p); + ++i; + } + return points; +} + + +// a = big radius, b = small radius + +template <typename Kernel> +std::vector<typename Kernel::Point_d> +generate_points_on_klein_bottle_variant_5D( + std::size_t num_points, double a, double b, bool uniform = false) { + typedef typename Kernel::Point_d Point; + typedef typename Kernel::FT FT; + Kernel k; + CGAL::Random rng; + + // if uniform + std::size_t num_lines = (std::size_t)sqrt(num_points); + + std::vector<Point> points; + points.reserve(num_points); + for (std::size_t i = 0; i < num_points;) { + FT u, v; + if (uniform) { + std::size_t k1 = i / num_lines; + std::size_t k2 = i % num_lines; + u = 6.2832 * k1 / num_lines; + v = 6.2832 * k2 / num_lines; + } else { + u = rng.get_double(0, 6.2832); + v = rng.get_double(0, 6.2832); + } + FT x1 = (a + b * cos(v)) * cos(u); + FT x2 = (a + b * cos(v)) * sin(u); + FT x3 = b * sin(v) * cos(u / 2); + FT x4 = b * sin(v) * sin(u / 2); + FT x5 = x1 + x2 + x3 + x4; + + Point p = construct_point(k, x1, x2, x3, x4, x5); + points.push_back(p); + ++i; + } + return points; +} + +} // namespace Gudhi + +#endif // RANDOM_POINT_GENERATORS_H_ diff --git a/src/common/include/gudhi/reader_utils.h b/src/common/include/gudhi/reader_utils.h index 899f9df6..97a87edd 100644 --- a/src/common/include/gudhi/reader_utils.h +++ b/src/common/include/gudhi/reader_utils.h @@ -2,9 +2,9 @@ * (Geometric Understanding in Higher Dimensions) is a generic C++ * library for computational topology. * - * Author(s): Clément Maria + * Author(s): Clement Maria, Pawel Dlotko * - * Copyright (C) 2014 INRIA Sophia Antipolis-Méditerranée (France) + * Copyright (C) 2014 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -30,18 +30,25 @@ #include <iostream> #include <fstream> #include <map> -#include <limits> // for numeric_limits<> +#include <limits> // for numeric_limits #include <string> #include <vector> +#include <utility> // for pair + +// Keep this file tag for Doxygen to parse the code, otherwise, functions are not documented. +// It is required for global functions and variables. + +/** @file + * @brief This file includes common file reader for GUDHI + */ /** - * \brief Read a set of points to turn it - * into a vector< vector<double> > by filling points + * @brief Read a set of points to turn it into a vector< vector<double> > by filling points. * - * File format: 1 point per line - * X11 X12 ... X1d - * X21 X22 ... X2d - * etc + * File format: 1 point per line<br> + * X11 X12 ... X1d<br> + * X21 X22 ... X2d<br> + * etc<br> */ inline void read_points(std::string file_name, std::vector< std::vector< double > > & points) { std::ifstream in_file(file_name.c_str(), std::ios::in); @@ -66,23 +73,29 @@ inline void read_points(std::string file_name, std::vector< std::vector< double } /** - * \brief Read a graph from a file. + * @brief Read a graph from a file. + * + * \tparam Graph_t Type for the return graph. Must be constructible from iterators on pairs of Vertex_handle + * \tparam Filtration_value Type for the value of the read filtration + * \tparam Vertex_handle Type for the value of the read vertices * - * File format: 1 simplex per line - * Dim1 X11 X12 ... X1d Fil1 - * Dim2 X21 X22 ... X2d Fil2 - * etc + * File format: 1 simplex per line<br> + * Dim1 X11 X12 ... X1d Fil1<br> + * Dim2 X21 X22 ... X2d Fil2<br> + * etc<br> * * The vertices must be labeled from 0 to n-1. * Every simplex must appear exactly once. * Simplices of dimension more than 1 are ignored. */ -inline Graph_t read_graph(std::string file_name) { +template< typename Graph_t, typename Filtration_value, typename Vertex_handle > +Graph_t read_graph(std::string file_name) { std::ifstream in_(file_name.c_str(), std::ios::in); if (!in_.is_open()) { std::cerr << "Unable to open file " << file_name << std::endl; } + typedef std::pair< Vertex_handle, Vertex_handle > Edge_t; std::vector< Edge_t > edges; std::vector< Filtration_value > edges_fil; std::map< Vertex_handle, Filtration_value > vertices; @@ -130,7 +143,7 @@ inline Graph_t read_graph(std::string file_name) { Graph_t skel_graph(edges.begin(), edges.end(), edges_fil.begin(), vertices.size()); auto vertex_prop = boost::get(vertex_filtration_t(), skel_graph); - boost::graph_traits<Graph_t>::vertex_iterator vi, vi_end; + typename boost::graph_traits<Graph_t>::vertex_iterator vi, vi_end; auto v_it = vertices.begin(); for (std::tie(vi, vi_end) = boost::vertices(skel_graph); vi != vi_end; ++vi, ++v_it) { boost::put(vertex_prop, *vi, v_it->second); @@ -140,12 +153,12 @@ inline Graph_t read_graph(std::string file_name) { } /** - * \brief Read a face from a file. + * @brief Read a face from a file. * - * File format: 1 simplex per line - * Dim1 X11 X12 ... X1d Fil1 - * Dim2 X21 X22 ... X2d Fil2 - * etc + * File format: 1 simplex per line<br> + * Dim1 X11 X12 ... X1d Fil1<br> + * Dim2 X21 X22 ... X2d Fil2<br> + * etc<br> * * The vertices must be labeled from 0 to n-1. * Every simplex must appear exactly once. @@ -166,18 +179,16 @@ bool read_simplex(std::istream & in_, std::vector< Vertex_handle > & simplex, Fi } /** - * \brief Read a hasse simplex from a file. - * - * File format: 1 simplex per line - * Dim1 k11 k12 ... k1Dim1 Fil1 - * Dim2 k21 k22 ... k2Dim2 Fil2 - * etc - * - * The key of a simplex is its position in the filtration order - * and also the number of its row in the file. - * Dimi ki1 ki2 ... kiDimi Fili means that the ith simplex in the - * filtration has dimension Dimi, filtration value fil1 and simplices with - * key ki1 ... kiDimi in its boundary.*/ + * @brief Read a hasse simplex from a file. + * + * File format: 1 simplex per line<br> + * Dim1 k11 k12 ... k1Dim1 Fil1<br> + * Dim2 k21 k22 ... k2Dim2 Fil2<br> + * etc<br> + * + * The key of a simplex is its position in the filtration order and also the number of its row in the file. + * Dimi ki1 ki2 ... kiDimi Fili means that the ith simplex in the filtration has dimension Dimi, filtration value + * fil1 and simplices with key ki1 ... kiDimi in its boundary.*/ template< typename Simplex_key, typename Filtration_value > bool read_hasse_simplex(std::istream & in_, std::vector< Simplex_key > & boundary, Filtration_value & fil) { int dim; @@ -195,4 +206,93 @@ bool read_hasse_simplex(std::istream & in_, std::vector< Simplex_key > & boundar return true; } +/** + * @brief Read a lower triangular distance matrix from a csv file. We assume that the .csv store the whole + * (square) matrix. + * + * @author Pawel Dlotko + * + * Square matrix file format:<br> + * 0;D12;...;D1j<br> + * D21;0;...;D2j<br> + * ...<br> + * Dj1;Dj2;...;0<br> + * + * lower matrix file format:<br> + * 0<br> + * D21;<br> + * D31;D32;<br> + * ...<br> + * Dj1;Dj2;...;Dj(j-1);<br> + * + **/ +template< typename Filtration_value > +std::vector< std::vector< Filtration_value > > read_lower_triangular_matrix_from_csv_file(const std::string& filename, + const char separator = ';') { +#ifdef DEBUG_TRACES + std::cout << "Using procedure read_lower_triangular_matrix_from_csv_file \n"; +#endif // DEBUG_TRACES + std::vector< std::vector< Filtration_value > > result; + std::ifstream in; + in.open(filename.c_str()); + if (!in.is_open()) { + return result; + } + + std::string line; + + // the first line is emtpy, so we ignore it: + std::getline(in, line); + std::vector< Filtration_value > values_in_this_line; + result.push_back(values_in_this_line); + + int number_of_line = 0; + + // first, read the file line by line to a string: + while (std::getline(in, line)) { + // if line is empty, break + if (line.size() == 0) + break; + + // if the last element of a string is comma: + if (line[ line.size() - 1 ] == separator) { + // then shrink the string by one + line.pop_back(); + } + + // replace all commas with spaces + std::replace(line.begin(), line.end(), separator, ' '); + + // put the new line to a stream + std::istringstream iss(line); + // and now read the doubles. + + int number_of_entry = 0; + std::vector< Filtration_value > values_in_this_line; + while (iss.good()) { + double entry; + iss >> entry; + if (number_of_entry <= number_of_line) { + values_in_this_line.push_back(entry); + } + ++number_of_entry; + } + if (!values_in_this_line.empty())result.push_back(values_in_this_line); + ++number_of_line; + } + in.close(); + +#ifdef DEBUG_TRACES + std::cerr << "Here is the matrix we read : \n"; + for (size_t i = 0; i != result.size(); ++i) { + for (size_t j = 0; j != result[i].size(); ++j) { + std::cerr << result[i][j] << " "; + } + std::cerr << std::endl; + } +#endif // DEBUG_TRACES + + return result; +} // read_lower_triangular_matrix_from_csv_file + #endif // READER_UTILS_H_ |