diff options
Diffstat (limited to 'include/gudhi')
100 files changed, 1505 insertions, 283 deletions
diff --git a/include/gudhi/Active_witness/Active_witness.h b/include/gudhi/Active_witness/Active_witness.h index d41a6811..8cb8662b 100644 --- a/include/gudhi/Active_witness/Active_witness.h +++ b/include/gudhi/Active_witness/Active_witness.h @@ -4,7 +4,7 @@ * * Author(s): Siargey Kachanovich * - * Copyright (C) 2016 INRIA (France) + * 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 diff --git a/include/gudhi/Active_witness/Active_witness_iterator.h b/include/gudhi/Active_witness/Active_witness_iterator.h index 0a05173a..10d2ec52 100644 --- a/include/gudhi/Active_witness/Active_witness_iterator.h +++ b/include/gudhi/Active_witness/Active_witness_iterator.h @@ -4,7 +4,7 @@ * * Author(s): Siargey Kachanovich * - * Copyright (C) 2016 INRIA (France) + * 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 diff --git a/include/gudhi/Alpha_complex.h b/include/gudhi/Alpha_complex.h index 63c6675c..4c07eddb 100644 --- a/include/gudhi/Alpha_complex.h +++ b/include/gudhi/Alpha_complex.h @@ -4,7 +4,7 @@ * * Author(s): Vincent Rouvreau * - * Copyright (C) 2015 INRIA + * Copyright (C) 2015 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 @@ -34,6 +34,7 @@ #include <CGAL/Epick_d.h> #include <CGAL/Spatial_sort_traits_adapter_d.h> #include <CGAL/property_map.h> // for CGAL::Identity_property_map +#include <CGAL/NT_converter.h> #include <iostream> #include <vector> @@ -323,8 +324,9 @@ class Alpha_complex { if (f_simplex_dim > 0) { // squared_radius function initialization Squared_Radius squared_radius = kernel_.compute_squared_radius_d_object(); + CGAL::NT_converter<typename Geom_traits::FT, Filtration_value> cv; - alpha_complex_filtration = squared_radius(pointVector.begin(), pointVector.end()); + alpha_complex_filtration = cv(squared_radius(pointVector.begin(), pointVector.end())); } complex.assign_filtration(f_simplex, alpha_complex_filtration); #ifdef DEBUG_TRACES diff --git a/include/gudhi/Bitmap_cubical_complex.h b/include/gudhi/Bitmap_cubical_complex.h index 969daba6..cc19b8b5 100644 --- a/include/gudhi/Bitmap_cubical_complex.h +++ b/include/gudhi/Bitmap_cubical_complex.h @@ -4,7 +4,7 @@ * * Author(s): Pawel Dlotko * - * Copyright (C) 2015 INRIA Sophia-Saclay (France) + * Copyright (C) 2015 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 @@ -383,7 +383,7 @@ class Bitmap_cubical_complex : public T { std::vector<std::size_t> bdry = this->get_boundary_of_a_cell(sh); if (globalDbg) { std::cerr << "std::pair<Simplex_handle, Simplex_handle> endpoints( Simplex_handle sh )\n"; - std::cerr << "bdry.size() : " << bdry.size() << std::endl; + std::cerr << "bdry.size() : " << bdry.size() << "\n"; } // this method returns two first elements from the boundary of sh. if (bdry.size() < 2) diff --git a/include/gudhi/Bitmap_cubical_complex/counter.h b/include/gudhi/Bitmap_cubical_complex/counter.h index 705b68a0..f82d4cc3 100644 --- a/include/gudhi/Bitmap_cubical_complex/counter.h +++ b/include/gudhi/Bitmap_cubical_complex/counter.h @@ -4,7 +4,7 @@ * * Author(s): Pawel Dlotko * - * Copyright (C) 2015 INRIA Sophia-Saclay (France) + * Copyright (C) 2015 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 diff --git a/include/gudhi/Bitmap_cubical_complex_base.h b/include/gudhi/Bitmap_cubical_complex_base.h index bf257be1..9b74e267 100644 --- a/include/gudhi/Bitmap_cubical_complex_base.h +++ b/include/gudhi/Bitmap_cubical_complex_base.h @@ -4,7 +4,7 @@ * * Author(s): Pawel Dlotko * - * Copyright (C) 2015 INRIA Sophia-Saclay (France) + * Copyright (C) 2015 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 @@ -492,7 +492,7 @@ class Bitmap_cubical_complex_base { this->multipliers.push_back(multiplier); multiplier *= 2 * sizes[i] + 1; } - this->data = std::vector<T>(multiplier, std::numeric_limits<T>::max()); + this->data = std::vector<T>(multiplier, std::numeric_limits<T>::infinity()); this->total_number_of_cells = multiplier; } @@ -562,7 +562,7 @@ void Bitmap_cubical_complex_base<T>::put_data_to_bins(T diameter_of_bin) { template <typename T> std::pair<T, T> Bitmap_cubical_complex_base<T>::min_max_filtration() { - std::pair<T, T> min_max(std::numeric_limits<T>::max(), std::numeric_limits<T>::min()); + std::pair<T, T> min_max(std::numeric_limits<T>::infinity(), -std::numeric_limits<T>::infinity()); for (std::size_t i = 0; i != this->data.size(); ++i) { if (this->data[i] < min_max.first) min_max.first = this->data[i]; if (this->data[i] > min_max.second) min_max.second = this->data[i]; diff --git a/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h b/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h index 4a0d1c74..8c35f590 100644 --- a/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h +++ b/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h @@ -4,7 +4,7 @@ * * Author(s): Pawel Dlotko * - * Copyright (C) 2015 INRIA Sophia-Saclay (France) + * Copyright (C) 2015 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 @@ -177,7 +177,7 @@ class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_c } } // std::reverse( this->sizes.begin() , this->sizes.end() ); - this->data = std::vector<T>(multiplier, std::numeric_limits<T>::max()); + this->data = std::vector<T>(multiplier, std::numeric_limits<T>::infinity()); this->total_number_of_cells = multiplier; } Bitmap_cubical_complex_periodic_boundary_conditions_base(const std::vector<unsigned>& sizes); diff --git a/include/gudhi/Bottleneck.h b/include/gudhi/Bottleneck.h index 7aee07bb..b0fc3949 100644 --- a/include/gudhi/Bottleneck.h +++ b/include/gudhi/Bottleneck.h @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA + * Copyright (C) 2015 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,6 +30,7 @@ #include <limits> // for numeric_limits #include <cmath> +#include <cfloat> // FLT_EVAL_METHOD namespace Gudhi { @@ -43,6 +44,13 @@ double bottleneck_distance_approx(Persistence_graph& g, double e) { Graph_matching biggest_unperfect(g); while (b_upper_bound - b_lower_bound > 2 * e) { double step = b_lower_bound + (b_upper_bound - b_lower_bound) / alpha; +#if !defined FLT_EVAL_METHOD || FLT_EVAL_METHOD < 0 || FLT_EVAL_METHOD > 1 + // On platforms where double computation is done with excess precision, + // we force it to its true precision so the following test is reliable. + volatile double drop_excess_precision = step; + step = drop_excess_precision; + // Alternative: step = CGAL::IA_force_to_double(step); +#endif if (step <= b_lower_bound || step >= b_upper_bound) // Avoid precision problem break; m.set_r(step); diff --git a/include/gudhi/Cech_complex.h b/include/gudhi/Cech_complex.h new file mode 100644 index 00000000..f9b8a269 --- /dev/null +++ b/include/gudhi/Cech_complex.h @@ -0,0 +1,130 @@ +/* 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): Vincent Rouvreau + * + * Copyright (C) 2018 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 CECH_COMPLEX_H_ +#define CECH_COMPLEX_H_ + +#include <gudhi/distance_functions.h> // for Gudhi::Minimal_enclosing_ball_radius +#include <gudhi/graph_simplicial_complex.h> // for Gudhi::Proximity_graph +#include <gudhi/Debug_utils.h> // for GUDHI_CHECK +#include <gudhi/Cech_complex_blocker.h> // for Gudhi::cech_complex::Cech_blocker + +#include <iostream> +#include <stdexcept> // for exception management +#include <vector> + +namespace Gudhi { + +namespace cech_complex { + +/** + * \class Cech_complex + * \brief Cech complex data structure. + * + * \ingroup cech_complex + * + * \details + * The data structure is a proximity graph, containing edges when the edge length is less or equal + * to a given max_radius. Edge length is computed from `Gudhi::Minimal_enclosing_ball_radius` distance function. + * + * \tparam SimplicialComplexForProximityGraph furnishes `Vertex_handle` and `Filtration_value` type definition required + * by `Gudhi::Proximity_graph`. + * + * \tparam ForwardPointRange must be a range for which `std::begin()` and `std::end()` methods return input + * iterators on a point. `std::begin()` and `std::end()` methods are also required for a point. + */ +template <typename SimplicialComplexForProximityGraph, typename ForwardPointRange> +class Cech_complex { + private: + // Required by compute_proximity_graph + using Vertex_handle = typename SimplicialComplexForProximityGraph::Vertex_handle; + using Filtration_value = typename SimplicialComplexForProximityGraph::Filtration_value; + using Proximity_graph = Gudhi::Proximity_graph<SimplicialComplexForProximityGraph>; + + // Retrieve Coordinate type from ForwardPointRange + using Point_from_range_iterator = typename boost::range_const_iterator<ForwardPointRange>::type; + using Point_from_range = typename std::iterator_traits<Point_from_range_iterator>::value_type; + using Coordinate_iterator = typename boost::range_const_iterator<Point_from_range>::type; + using Coordinate = typename std::iterator_traits<Coordinate_iterator>::value_type; + + public: + // Point and Point_cloud type definition + using Point = std::vector<Coordinate>; + using Point_cloud = std::vector<Point>; + + public: + /** \brief Cech_complex constructor from a list of points. + * + * @param[in] points Range of points. + * @param[in] max_radius Maximal radius value. + * + * \tparam ForwardPointRange must be a range of Point. Point must be a range of <b>copyable</b> Cartesian coordinates. + * + */ + Cech_complex(const ForwardPointRange& points, Filtration_value max_radius) : max_radius_(max_radius) { + // Point cloud deep copy + point_cloud_.reserve(boost::size(points)); + for (auto&& point : points) point_cloud_.emplace_back(std::begin(point), std::end(point)); + + cech_skeleton_graph_ = Gudhi::compute_proximity_graph<SimplicialComplexForProximityGraph>( + point_cloud_, max_radius_, Gudhi::Minimal_enclosing_ball_radius()); + } + + /** \brief Initializes the simplicial complex from the proximity graph and expands it until a given maximal + * dimension, using the Cech blocker oracle. + * + * @param[in] complex SimplicialComplexForCech to be created. + * @param[in] dim_max graph expansion until this given maximal dimension. + * @exception std::invalid_argument In debug mode, if `complex.num_vertices()` does not return 0. + * + */ + template <typename SimplicialComplexForCechComplex> + void create_complex(SimplicialComplexForCechComplex& complex, int dim_max) { + GUDHI_CHECK(complex.num_vertices() == 0, + std::invalid_argument("Cech_complex::create_complex - simplicial complex is not empty")); + + // insert the proximity graph in the simplicial complex + complex.insert_graph(cech_skeleton_graph_); + // expand the graph until dimension dim_max + complex.expansion_with_blockers(dim_max, + Cech_blocker<SimplicialComplexForCechComplex, Cech_complex>(&complex, this)); + } + + /** @return max_radius value given at construction. */ + Filtration_value max_radius() const { return max_radius_; } + + /** @param[in] vertex Point position in the range. + * @return The point. + */ + const Point& get_point(Vertex_handle vertex) const { return point_cloud_[vertex]; } + + private: + Proximity_graph cech_skeleton_graph_; + Filtration_value max_radius_; + Point_cloud point_cloud_; +}; + +} // namespace cech_complex + +} // namespace Gudhi + +#endif // CECH_COMPLEX_H_ diff --git a/include/gudhi/Cech_complex_blocker.h b/include/gudhi/Cech_complex_blocker.h new file mode 100644 index 00000000..b0d347b1 --- /dev/null +++ b/include/gudhi/Cech_complex_blocker.h @@ -0,0 +1,91 @@ +/* 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): Vincent Rouvreau + * + * Copyright (C) 2018 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 CECH_COMPLEX_BLOCKER_H_ +#define CECH_COMPLEX_BLOCKER_H_ + +#include <gudhi/distance_functions.h> // for Gudhi::Minimal_enclosing_ball_radius + +#include <iostream> +#include <vector> +#include <cmath> // for std::sqrt + +namespace Gudhi { + +namespace cech_complex { + +/** \internal + * \class Cech_blocker + * \brief Čech complex blocker. + * + * \ingroup cech_complex + * + * \details + * Čech blocker is an oracle constructed from a Cech_complex and a simplicial complex. + * + * \tparam SimplicialComplexForProximityGraph furnishes `Simplex_handle` and `Filtration_value` type definition, + * `simplex_vertex_range(Simplex_handle sh)`and `assign_filtration(Simplex_handle sh, Filtration_value filt)` methods. + * + * \tparam Chech_complex is required by the blocker. + */ +template <typename SimplicialComplexForCech, typename Cech_complex> +class Cech_blocker { + private: + using Point_cloud = typename Cech_complex::Point_cloud; + + using Simplex_handle = typename SimplicialComplexForCech::Simplex_handle; + using Filtration_value = typename SimplicialComplexForCech::Filtration_value; + + public: + /** \internal \brief Čech complex blocker operator() - the oracle - assigns the filtration value from the simplex + * radius and returns if the simplex expansion must be blocked. + * \param[in] sh The Simplex_handle. + * \return true if the simplex radius is greater than the Cech_complex max_radius*/ + bool operator()(Simplex_handle sh) { + Point_cloud points; + for (auto vertex : sc_ptr_->simplex_vertex_range(sh)) { + points.push_back(cc_ptr_->get_point(vertex)); +#ifdef DEBUG_TRACES + std::cout << "#(" << vertex << ")#"; +#endif // DEBUG_TRACES + } + Filtration_value radius = Gudhi::Minimal_enclosing_ball_radius()(points); +#ifdef DEBUG_TRACES + if (radius > cc_ptr_->max_radius()) std::cout << "radius > max_radius => expansion is blocked\n"; +#endif // DEBUG_TRACES + sc_ptr_->assign_filtration(sh, radius); + return (radius > cc_ptr_->max_radius()); + } + + /** \internal \brief Čech complex blocker constructor. */ + Cech_blocker(SimplicialComplexForCech* sc_ptr, Cech_complex* cc_ptr) : sc_ptr_(sc_ptr), cc_ptr_(cc_ptr) {} + + private: + SimplicialComplexForCech* sc_ptr_; + Cech_complex* cc_ptr_; +}; + +} // namespace cech_complex + +} // namespace Gudhi + +#endif // CECH_COMPLEX_BLOCKER_H_ diff --git a/include/gudhi/Clock.h b/include/gudhi/Clock.h index b83de2f5..cdf18cb2 100644 --- a/include/gudhi/Clock.h +++ b/include/gudhi/Clock.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 diff --git a/include/gudhi/Contraction/Edge_profile.h b/include/gudhi/Contraction/Edge_profile.h index e4910b27..30b1b80a 100644 --- a/include/gudhi/Contraction/Edge_profile.h +++ b/include/gudhi/Contraction/Edge_profile.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 diff --git a/include/gudhi/Contraction/policies/Contraction_visitor.h b/include/gudhi/Contraction/policies/Contraction_visitor.h index 7ee05aad..fa02308b 100644 --- a/include/gudhi/Contraction/policies/Contraction_visitor.h +++ b/include/gudhi/Contraction/policies/Contraction_visitor.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 diff --git a/include/gudhi/Contraction/policies/Cost_policy.h b/include/gudhi/Contraction/policies/Cost_policy.h index f4d343ec..04ce36b6 100644 --- a/include/gudhi/Contraction/policies/Cost_policy.h +++ b/include/gudhi/Contraction/policies/Cost_policy.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 diff --git a/include/gudhi/Contraction/policies/Dummy_valid_contraction.h b/include/gudhi/Contraction/policies/Dummy_valid_contraction.h index 5d329496..a5567454 100644 --- a/include/gudhi/Contraction/policies/Dummy_valid_contraction.h +++ b/include/gudhi/Contraction/policies/Dummy_valid_contraction.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 diff --git a/include/gudhi/Contraction/policies/Edge_length_cost.h b/include/gudhi/Contraction/policies/Edge_length_cost.h index dac2d448..1b7a825b 100644 --- a/include/gudhi/Contraction/policies/Edge_length_cost.h +++ b/include/gudhi/Contraction/policies/Edge_length_cost.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 diff --git a/include/gudhi/Contraction/policies/First_vertex_placement.h b/include/gudhi/Contraction/policies/First_vertex_placement.h index 1f68db0d..0b9f8775 100644 --- a/include/gudhi/Contraction/policies/First_vertex_placement.h +++ b/include/gudhi/Contraction/policies/First_vertex_placement.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 diff --git a/include/gudhi/Contraction/policies/Link_condition_valid_contraction.h b/include/gudhi/Contraction/policies/Link_condition_valid_contraction.h index 250bba27..8c869830 100644 --- a/include/gudhi/Contraction/policies/Link_condition_valid_contraction.h +++ b/include/gudhi/Contraction/policies/Link_condition_valid_contraction.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 diff --git a/include/gudhi/Contraction/policies/Middle_placement.h b/include/gudhi/Contraction/policies/Middle_placement.h index 4b59f1b5..0ba23a35 100644 --- a/include/gudhi/Contraction/policies/Middle_placement.h +++ b/include/gudhi/Contraction/policies/Middle_placement.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 diff --git a/include/gudhi/Contraction/policies/Placement_policy.h b/include/gudhi/Contraction/policies/Placement_policy.h index 34ffa49f..19509fad 100644 --- a/include/gudhi/Contraction/policies/Placement_policy.h +++ b/include/gudhi/Contraction/policies/Placement_policy.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 diff --git a/include/gudhi/Contraction/policies/Valid_contraction_policy.h b/include/gudhi/Contraction/policies/Valid_contraction_policy.h index 78d61173..8a91f0b5 100644 --- a/include/gudhi/Contraction/policies/Valid_contraction_policy.h +++ b/include/gudhi/Contraction/policies/Valid_contraction_policy.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 diff --git a/include/gudhi/Debug_utils.h b/include/gudhi/Debug_utils.h index 90d3cf47..3f5cb04f 100644 --- a/include/gudhi/Debug_utils.h +++ b/include/gudhi/Debug_utils.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/Edge_contraction.h b/include/gudhi/Edge_contraction.h index cf9a2c27..fcd06996 100644 --- a/include/gudhi/Edge_contraction.h +++ b/include/gudhi/Edge_contraction.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 diff --git a/include/gudhi/Euclidean_strong_witness_complex.h b/include/gudhi/Euclidean_strong_witness_complex.h index 4f3cef4f..ea97cd3f 100644 --- a/include/gudhi/Euclidean_strong_witness_complex.h +++ b/include/gudhi/Euclidean_strong_witness_complex.h @@ -4,7 +4,7 @@ * * Author(s): Siargey Kachanovich * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 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 diff --git a/include/gudhi/Euclidean_witness_complex.h b/include/gudhi/Euclidean_witness_complex.h index ff8bb139..1dacefa5 100644 --- a/include/gudhi/Euclidean_witness_complex.h +++ b/include/gudhi/Euclidean_witness_complex.h @@ -4,7 +4,7 @@ * * Author(s): Siargey Kachanovich * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 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 diff --git a/include/gudhi/GIC.h b/include/gudhi/GIC.h index 40ff7a4a..7aa95210 100644 --- a/include/gudhi/GIC.h +++ b/include/gudhi/GIC.h @@ -4,7 +4,7 @@ * * Author: Mathieu Carriere * - * Copyright (C) 2017 INRIA + * Copyright (C) 2017 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 @@ -23,6 +23,11 @@ #ifndef GIC_H_ #define GIC_H_ +#ifdef GUDHI_USE_TBB +#include <tbb/parallel_for.h> +#include <tbb/mutex.h> +#endif + #include <gudhi/Debug_utils.h> #include <gudhi/graph_simplicial_complex.h> #include <gudhi/reader_utils.h> @@ -99,9 +104,8 @@ class Cover_complex { int data_dimension; // dimension of input data. int n; // number of points. - std::map<int, double> func; // function used to compute the output simplicial complex. - std::map<int, double> - func_color; // function used to compute the colors of the nodes of the output simplicial complex. + std::vector<double> func; // function used to compute the output simplicial complex. + std::vector<double> func_color; // function used to compute the colors of the nodes of the output simplicial complex. bool functional_cover = false; // whether we use a cover with preimages of a function or not. Graph one_skeleton_OFF; // one-skeleton given by the input OFF file (if it exists). @@ -114,8 +118,8 @@ class Cover_complex { Persistence_diagram PD; std::vector<double> distribution; - std::map<int, std::vector<int> > - cover; // function associating to each data point its vectors of cover elements to which it belongs. + std::vector<std::vector<int> > + cover; // function associating to each data point the vector of cover elements to which it belongs. std::map<int, std::vector<int> > cover_back; // inverse of cover, in order to get the data points associated to a specific cover element. std::map<int, double> cover_std; // standard function (induced by func) used to compute the extended persistence @@ -138,28 +142,26 @@ class Cover_complex { std::string point_cloud_name; std::string color_name; - // Point comparator - struct Less { - Less(std::map<int, double> func) { Fct = func; } - std::map<int, double> Fct; - bool operator()(int a, int b) { - if (Fct[a] == Fct[b]) - return a < b; - else - return Fct[a] < Fct[b]; - } - }; - // Remove all edges of a graph. void remove_edges(Graph& G) { boost::graph_traits<Graph>::edge_iterator ei, ei_end; for (boost::tie(ei, ei_end) = boost::edges(G); ei != ei_end; ++ei) boost::remove_edge(*ei, G); } + // Thread local is not available on XCode version < V.8 + // If not available, random engine is a class member. +#ifndef GUDHI_CAN_USE_CXX11_THREAD_LOCAL + std::default_random_engine re; +#endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL + // Find random number in [0,1]. double GetUniform() { + // Thread local is not available on XCode version < V.8 + // If available, random engine is defined for each thread. +#ifdef GUDHI_CAN_USE_CXX11_THREAD_LOCAL thread_local std::default_random_engine re; - thread_local std::uniform_real_distribution<double> Dist(0, 1); +#endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL + std::uniform_real_distribution<double> Dist(0, 1); return Dist(re); } @@ -276,6 +278,7 @@ class Cover_complex { point_cloud.emplace_back(point.begin(), point.begin() + data_dimension); boost::add_vertex(one_skeleton_OFF); vertices.push_back(boost::add_vertex(one_skeleton)); + cover.emplace_back(); i++; } } @@ -422,7 +425,6 @@ class Cover_complex { 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<int> samples(m); double delta = 0; if (verbose) std::cout << n << " points in R^" << data_dimension << std::endl; @@ -430,17 +432,36 @@ class Cover_complex { 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); + // This cannot be parallelized if thread_local is not defined + // thread_local is not defined for XCode < v.8 + #if defined(GUDHI_USE_TBB) && defined(GUDHI_CAN_USE_CXX11_THREAD_LOCAL) + tbb::mutex deltamutex; + tbb::parallel_for(0, N, [&](int i){ + std::vector<int> samples(m); + 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); + } + deltamutex.lock(); + delta += hausdorff_dist / N; + deltamutex.unlock(); + }); + #else + for (int i = 0; i < N; i++) { + std::vector<int> samples(m); + 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; } - delta += hausdorff_dist / N; - } + #endif if (verbose) std::cout << "delta = " << delta << std::endl; set_graph_from_rips(delta, distance); @@ -465,7 +486,7 @@ class Cover_complex { while (std::getline(input, line)) { std::stringstream stream(line); stream >> f; - func.emplace(i, f); + func.push_back(f); i++; } functional_cover = true; @@ -479,7 +500,7 @@ class Cover_complex { * */ void set_function_from_coordinate(int k) { - for (int i = 0; i < n; i++) func.emplace(i, point_cloud[i][k]); + for (int i = 0; i < n; i++) func.push_back(point_cloud[i][k]); functional_cover = true; cover_name = "coordinate " + std::to_string(k); } @@ -492,7 +513,7 @@ class Cover_complex { */ template <class InputRange> void set_function_from_range(InputRange const& function) { - for (int i = 0; i < n; i++) func.emplace(i, function[i]); + for (int i = 0; i < n; i++) func.push_back(function[i]); functional_cover = true; } @@ -654,7 +675,7 @@ class Cover_complex { // Sort points according to function values std::vector<int> points(n); for (int i = 0; i < n; i++) points[i] = i; - std::sort(points.begin(), points.end(), Less(this->func)); + std::sort(points.begin(), points.end(), [=](const int & p1, const int & p2){return (this->func[p1] < this->func[p2]);}); int id = 0; int pos = 0; @@ -710,37 +731,74 @@ class Cover_complex { funcstd[i] = 0.5 * (u + v); } - if (verbose) std::cout << "Computing connected components..." << std::endl; - // #pragma omp parallel for - for (int i = 0; i < res; i++) { - // Compute connected components - Graph G = one_skeleton.create_subgraph(); - int num = preimages[i].size(); - std::vector<int> component(num); - for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G); - boost::connected_components(G, &component[0]); - int max = 0; - - // For each point in preimage - for (int j = 0; j < num; j++) { - // Update number of components in preimage - if (component[j] > max) max = component[j]; - - // Identify component with Cantor polynomial N^2 -> N - int identifier = (std::pow(i + component[j], 2) + 3 * i + component[j]) / 2; - - // Update covers - cover[preimages[i][j]].push_back(identifier); - cover_back[identifier].push_back(preimages[i][j]); - cover_fct[identifier] = i; - cover_std[identifier] = funcstd[i]; - cover_color[identifier].second += func_color[preimages[i][j]]; - cover_color[identifier].first += 1; - } + #ifdef GUDHI_USE_TBB + if (verbose) std::cout << "Computing connected components (parallelized)..." << std::endl; + tbb::mutex covermutex, idmutex; + tbb::parallel_for(0, res, [&](int i){ + // Compute connected components + Graph G = one_skeleton.create_subgraph(); + int num = preimages[i].size(); + std::vector<int> component(num); + for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G); + boost::connected_components(G, &component[0]); + int max = 0; + + // For each point in preimage + for (int j = 0; j < num; j++) { + // Update number of components in preimage + if (component[j] > max) max = component[j]; + + // Identify component with Cantor polynomial N^2 -> N + int identifier = ((i + component[j])*(i + component[j]) + 3 * i + component[j]) / 2; + + // Update covers + covermutex.lock(); + cover[preimages[i][j]].push_back(identifier); + cover_back[identifier].push_back(preimages[i][j]); + cover_fct[identifier] = i; + cover_std[identifier] = funcstd[i]; + cover_color[identifier].second += func_color[preimages[i][j]]; + cover_color[identifier].first += 1; + covermutex.unlock(); + } - // Maximal dimension is total number of connected components - id += max + 1; - } + // Maximal dimension is total number of connected components + idmutex.lock(); + id += max + 1; + idmutex.unlock(); + }); + #else + if (verbose) std::cout << "Computing connected components..." << std::endl; + for (int i = 0; i < res; i++) { + // Compute connected components + Graph G = one_skeleton.create_subgraph(); + int num = preimages[i].size(); + std::vector<int> component(num); + for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G); + boost::connected_components(G, &component[0]); + int max = 0; + + // For each point in preimage + for (int j = 0; j < num; j++) { + // Update number of components in preimage + if (component[j] > max) max = component[j]; + + // Identify component with Cantor polynomial N^2 -> N + int identifier = (std::pow(i + component[j], 2) + 3 * i + component[j]) / 2; + + // Update covers + cover[preimages[i][j]].push_back(identifier); + cover_back[identifier].push_back(preimages[i][j]); + cover_fct[identifier] = i; + cover_std[identifier] = funcstd[i]; + cover_color[identifier].second += func_color[preimages[i][j]]; + cover_color[identifier].first += 1; + } + + // Maximal dimension is total number of connected components + id += max + 1; + } + #endif maximal_dim = id - 1; for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) @@ -803,24 +861,46 @@ class Cover_complex { for (int j = 0; j < n; j++) mindist[j] = std::numeric_limits<double>::max(); // Compute the geodesic distances to subsamples with Dijkstra - // #pragma omp parallel for - 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<double> dmap(n); - boost::dijkstra_shortest_paths( - one_skeleton, vertices[seed], - boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index))); - - for (int j = 0; j < n; j++) - if (mindist[j] > dmap[j]) { - mindist[j] = dmap[j]; - if (cover[j].size() == 0) - cover[j].push_back(i); - else - cover[j][0] = i; - } - } + #ifdef GUDHI_USE_TBB + if (verbose) std::cout << "Computing geodesic distances (parallelized)..." << std::endl; + tbb::mutex coverMutex; tbb::mutex mindistMutex; + tbb::parallel_for(0, m, [&](int i){ + int seed = voronoi_subsamples[i]; + std::vector<double> dmap(n); + boost::dijkstra_shortest_paths( + one_skeleton, vertices[seed], + boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index))); + + coverMutex.lock(); mindistMutex.lock(); + for (int j = 0; j < n; j++) + if (mindist[j] > dmap[j]) { + mindist[j] = dmap[j]; + if (cover[j].size() == 0) + cover[j].push_back(i); + else + cover[j][0] = i; + } + coverMutex.unlock(); mindistMutex.unlock(); + }); + #else + 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<double> dmap(n); + boost::dijkstra_shortest_paths( + one_skeleton, vertices[seed], + boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index))); + + for (int j = 0; j < n; j++) + if (mindist[j] > dmap[j]) { + mindist[j] = dmap[j]; + if (cover[j].size() == 0) + cover[j].push_back(i); + else + cover[j][0] = i; + } + } + #endif for (int i = 0; i < n; i++) { cover_back[cover[i][0]].push_back(i); @@ -860,7 +940,7 @@ class Cover_complex { while (std::getline(input, line)) { std::stringstream stream(line); stream >> f; - func_color.emplace(i, f); + func_color.push_back(f); i++; } color_name = color_file_name; @@ -873,7 +953,7 @@ class Cover_complex { * */ void set_color_from_coordinate(int k = 0) { - for (int i = 0; i < n; i++) func_color[i] = point_cloud[i][k]; + for (int i = 0; i < n; i++) func_color.push_back(point_cloud[i][k]); color_name = "coordinate "; color_name.append(std::to_string(k)); } @@ -885,7 +965,7 @@ class Cover_complex { * */ void set_color_from_vector(std::vector<double> color) { - for (unsigned int i = 0; i < color.size(); i++) func_color[i] = color[i]; + for (unsigned int i = 0; i < color.size(); i++) func_color.push_back(color[i]); } public: // Create a .dot file that can be compiled with neato to produce a .pdf file. @@ -1039,45 +1119,29 @@ class Cover_complex { minf = std::min(minf, it->second); } + // Build filtration for (auto const& simplex : simplices) { - // Add a simplex and a cone on it - std::vector<int> splx = simplex; - splx.push_back(-2); - st.insert_simplex_and_subfaces(splx); + std::vector<int> splx = simplex; splx.push_back(-2); + st.insert_simplex_and_subfaces(splx, -3); } - // Build filtration - for (auto simplex : st.complex_simplex_range()) { - double filta = std::numeric_limits<double>::lowest(); - double filts = filta; - bool ascending = true; - for (auto vertex : st.simplex_vertex_range(simplex)) { - if (vertex == -2) { - ascending = false; - continue; - } - filta = std::max(-2 + (cover_std[vertex] - minf) / (maxf - minf), filta); - filts = std::max(2 - (cover_std[vertex] - minf) / (maxf - minf), filts); - } - if (ascending) - st.assign_filtration(simplex, filta); - else - st.assign_filtration(simplex, filts); + for (std::map<int, double>::iterator it = cover_std.begin(); it != cover_std.end(); it++) { + int vertex = it->first; float val = it->second; + int vert[] = {vertex}; int edge[] = {vertex, -2}; + st.assign_filtration(st.find(vert), -2 + (val - minf)/(maxf - minf)); + st.assign_filtration(st.find(edge), 2 - (val - minf)/(maxf - minf)); } - int magic[] = {-2}; - st.assign_filtration(st.find(magic), -3); + st.make_filtration_non_decreasing(); // Compute PD - st.initialize_filtration(); - Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp> pcoh(st); - pcoh.init_coefficients(2); + Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp> pcoh(st); pcoh.init_coefficients(2); pcoh.compute_persistent_cohomology(); // Output PD int max_dim = st.dimension(); for (int i = 0; i < max_dim; i++) { std::vector<std::pair<double, double> > bars = pcoh.intervals_in_dimension(i); - int num_bars = bars.size(); + int num_bars = bars.size(); if(i == 0) num_bars -= 1; if(verbose) std::cout << num_bars << " interval(s) in dimension " << i << ":" << std::endl; for (int j = 0; j < num_bars; j++) { double birth = bars[j].first; @@ -1103,22 +1167,25 @@ class Cover_complex { * @param[in] N number of bootstrap iterations. * */ - template <typename SimplicialComplex> - void compute_distribution(int N = 100) { - if (distribution.size() >= N) { + void compute_distribution(unsigned int N = 100) { + unsigned int sz = distribution.size(); + if (sz >= N) { std::cout << "Already done!" << std::endl; } else { - for (int i = 0; i < N - distribution.size(); i++) { - Cover_complex Cboot; - Cboot.n = this->n; + for (unsigned int i = 0; i < N - sz; i++) { + if (verbose) std::cout << "Computing " << i << "th bootstrap, bottleneck distance = "; + + Cover_complex Cboot; Cboot.n = this->n; Cboot.data_dimension = this->data_dimension; Cboot.type = this->type; Cboot.functional_cover = true; + std::vector<int> boot(this->n); for (int j = 0; j < this->n; j++) { double u = GetUniform(); - int id = std::floor(u * (this->n)); - boot[j] = id; - Cboot.point_cloud[j] = this->point_cloud[id]; - Cboot.func.emplace(j, this->func[id]); + int id = std::floor(u * (this->n)); boot[j] = id; + Cboot.point_cloud.push_back(this->point_cloud[id]); Cboot.cover.emplace_back(); Cboot.func.push_back(this->func[id]); + boost::add_vertex(Cboot.one_skeleton_OFF); Cboot.vertices.push_back(boost::add_vertex(Cboot.one_skeleton)); } + Cboot.set_color_from_vector(Cboot.func); + for (int j = 0; j < n; j++) { std::vector<double> dist(n); for (int k = 0; k < n; k++) dist[k] = distances[boot[j]][boot[k]]; @@ -1131,8 +1198,9 @@ class Cover_complex { Cboot.set_cover_from_function(); Cboot.find_simplices(); Cboot.compute_PD(); - - distribution.push_back(Gudhi::persistence_diagram::bottleneck_distance(this->PD, Cboot.PD)); + double db = Gudhi::persistence_diagram::bottleneck_distance(this->PD, Cboot.PD); + if (verbose) std::cout << db << std::endl; + distribution.push_back(db); } std::sort(distribution.begin(), distribution.end()); @@ -1146,7 +1214,7 @@ class Cover_complex { * */ double compute_distance_from_confidence_level(double alpha) { - int N = distribution.size(); + unsigned int N = distribution.size(); return distribution[std::floor(alpha * N)]; } @@ -1157,9 +1225,11 @@ class Cover_complex { * */ double compute_confidence_level_from_distance(double d) { - int N = distribution.size(); - for (int i = 0; i < N; i++) - if (distribution[i] > d) return i * 1.0 / N; + unsigned int N = distribution.size(); + double level = 1; + for (unsigned int i = 0; i < N; i++) + if (distribution[i] > d){ level = i * 1.0 / N; break; } + return level; } public: @@ -1171,7 +1241,9 @@ class Cover_complex { double distancemin = -std::numeric_limits<double>::lowest(); int N = PD.size(); for (int i = 0; i < N; i++) distancemin = std::min(distancemin, 0.5 * (PD[i].second - PD[i].first)); - return 1 - compute_confidence_level_from_distance(distancemin); + double p_value = 1 - compute_confidence_level_from_distance(distancemin); + if (verbose) std::cout << "p value = " << p_value << std::endl; + return p_value; } // ******************************************************************************************************************* @@ -1206,8 +1278,7 @@ class Cover_complex { } if (type == "Nerve") { - for(auto& simplex : cover) - simplices.push_back(simplex.second); + for(int i = 0; i < n; i++) simplices.push_back(cover[i]); std::sort(simplices.begin(), simplices.end()); std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end()); simplices.resize(std::distance(simplices.begin(), it)); diff --git a/include/gudhi/Graph_matching.h b/include/gudhi/Graph_matching.h index f51e22e9..313e7d9c 100644 --- a/include/gudhi/Graph_matching.h +++ b/include/gudhi/Graph_matching.h @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA + * Copyright (C) 2015 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 diff --git a/include/gudhi/Hasse_complex.h b/include/gudhi/Hasse_complex.h index e67f7609..efcaea55 100644 --- a/include/gudhi/Hasse_complex.h +++ b/include/gudhi/Hasse_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 diff --git a/include/gudhi/Internal_point.h b/include/gudhi/Internal_point.h index 0b2d26fe..7f350f64 100644 --- a/include/gudhi/Internal_point.h +++ b/include/gudhi/Internal_point.h @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA + * Copyright (C) 2015 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 diff --git a/include/gudhi/Kd_tree_search.h b/include/gudhi/Kd_tree_search.h index 96bbeb36..ad1054e5 100644 --- a/include/gudhi/Kd_tree_search.h +++ b/include/gudhi/Kd_tree_search.h @@ -4,7 +4,7 @@ * * Author(s): Clement Jamin * - * Copyright (C) 2016 INRIA + * 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 diff --git a/include/gudhi/Miniball.COPYRIGHT b/include/gudhi/Miniball.COPYRIGHT new file mode 100644 index 00000000..dbe4c553 --- /dev/null +++ b/include/gudhi/Miniball.COPYRIGHT @@ -0,0 +1,4 @@ +The miniball software is available under the GNU General Public License (GPLv3 - https://www.gnu.org/copyleft/gpl.html). +If your intended use is not compliant with this license, please buy a commercial license (EUR 500 - https://people.inf.ethz.ch/gaertner/subdir/software/miniball/license.html). +You need a license if the software that you develop using Miniball V3.0 is not open source. + diff --git a/include/gudhi/Miniball.README b/include/gudhi/Miniball.README new file mode 100644 index 00000000..033d8953 --- /dev/null +++ b/include/gudhi/Miniball.README @@ -0,0 +1,26 @@ +https://people.inf.ethz.ch/gaertner/subdir/software/miniball.html + +Smallest Enclosing Balls of Points - Fast and Robust in C++. +(high-quality software for smallest enclosing balls of balls is available in the computational geometry algorithms library CGAL) + + +This is the miniball software (V3.0) for computing smallest enclosing balls of points in arbitrary dimensions. It consists of a C++ header file Miniball.hpp (around 500 lines of code) and two example programs miniball_example.cpp and miniball_example_containers.cpp that demonstrate the usage. The first example stores the coordinates of the input points in a two-dimensional array, the second example uses a list of vectors to show how generic containers can be used. + +Credits: Aditya Gupta and Alexandros Konstantinakis-Karmis have significantly contributed to this version of the software. + +Changes - https://people.inf.ethz.ch/gaertner/subdir/software/miniball/changes.txt - from previous versions. + +The theory - https://people.inf.ethz.ch/gaertner/subdir/texts/own_work/esa99_final.pdf - behind the miniball software (Proc. 7th Annual European Symposium on Algorithms (ESA), Lecture Notes in Computer Science 1643, Springer-Verlag, pp.325-338, 1999). + +Main Features: + + Very fast in low dimensions. 1 million points in 5-space are processed within 0.05 seconds on any recent machine. + + High numerical stability. Almost all input degeneracies (cospherical points, multiple points, points very close together) are routinely handled. + + Easily integrates into your code. You can freely choose the coordinate type of your points and the container to store the points. If you still need to adapt the code, the header is small and readable and contains documentation for all major methods. + + +Changes done for the GUDHI version of MiniBall: + - Add include guard + - Move Miniball namespace inside a new Gudhi namespace diff --git a/include/gudhi/Miniball.hpp b/include/gudhi/Miniball.hpp new file mode 100644 index 00000000..ce6cbb5b --- /dev/null +++ b/include/gudhi/Miniball.hpp @@ -0,0 +1,523 @@ +// Copright (C) 1999-2013, Bernd Gaertner +// $Rev: 3581 $ +// +// 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/>. +// +// Contact: +// -------- +// Bernd Gaertner +// Institute of Theoretical Computer Science +// ETH Zuerich +// CAB G31.1 +// CH-8092 Zuerich, Switzerland +// http://www.inf.ethz.ch/personal/gaertner + +#ifndef MINIBALL_HPP_ +#define MINIBALL_HPP_ + +#include <cassert> +#include <algorithm> +#include <list> +#include <ctime> +#include <limits> + +namespace Gudhi { + +namespace Miniball { + + // Global Functions + // ================ + template <typename NT> + inline NT mb_sqr (NT r) {return r*r;} + + // Functors + // ======== + + // functor to map a point iterator to the corresponding coordinate iterator; + // generic version for points whose coordinate containers have begin() + template < typename Pit_, typename Cit_ > + struct CoordAccessor { + typedef Pit_ Pit; + typedef Cit_ Cit; + inline Cit operator() (Pit it) const { return (*it).begin(); } + }; + + // partial specialization for points whose coordinate containers are arrays + template < typename Pit_, typename Cit_ > + struct CoordAccessor<Pit_, Cit_*> { + typedef Pit_ Pit; + typedef Cit_* Cit; + inline Cit operator() (Pit it) const { return *it; } + }; + + // Class Declaration + // ================= + + template <typename CoordAccessor> + class Miniball { + private: + // types + // The iterator type to go through the input points + typedef typename CoordAccessor::Pit Pit; + // The iterator type to go through the coordinates of a single point. + typedef typename CoordAccessor::Cit Cit; + // The coordinate type + typedef typename std::iterator_traits<Cit>::value_type NT; + // The iterator to go through the support points + typedef typename std::list<Pit>::iterator Sit; + + // data members... + const int d; // dimension + Pit points_begin; + Pit points_end; + CoordAccessor coord_accessor; + double time; + const NT nt0; // NT(0) + + //...for the algorithms + std::list<Pit> L; + Sit support_end; + int fsize; // number of forced points + int ssize; // number of support points + + // ...for the ball updates + NT* current_c; + NT current_sqr_r; + NT** c; + NT* sqr_r; + + // helper arrays + NT* q0; + NT* z; + NT* f; + NT** v; + NT** a; + + public: + // The iterator type to go through the support points + typedef typename std::list<Pit>::const_iterator SupportPointIterator; + + // PRE: [begin, end) is a nonempty range + // POST: computes the smallest enclosing ball of the points in the range + // [begin, end); the functor a maps a point iterator to an iterator + // through the d coordinates of the point + Miniball (int d_, Pit begin, Pit end, CoordAccessor ca = CoordAccessor()); + + // POST: returns a pointer to the first element of an array that holds + // the d coordinates of the center of the computed ball + const NT* center () const; + + // POST: returns the squared radius of the computed ball + NT squared_radius () const; + + // POST: returns the number of support points of the computed ball; + // the support points form a minimal set with the same smallest + // enclosing ball as the input set; in particular, the support + // points are on the boundary of the computed ball, and their + // number is at most d+1 + int nr_support_points () const; + + // POST: returns an iterator to the first support point + SupportPointIterator support_points_begin () const; + + // POST: returns a past-the-end iterator for the range of support points + SupportPointIterator support_points_end () const; + + // POST: returns the maximum excess of any input point w.r.t. the computed + // ball, divided by the squared radius of the computed ball. The + // excess of a point is the difference between its squared distance + // from the center and the squared radius; Ideally, the return value + // is 0. subopt is set to the absolute value of the most negative + // coefficient in the affine combination of the support points that + // yields the center. Ideally, this is a convex combination, and there + // is no negative coefficient in which case subopt is set to 0. + NT relative_error (NT& subopt) const; + + // POST: return true if the relative error is at most tol, and the + // suboptimality is 0; the default tolerance is 10 times the + // coordinate type's machine epsilon + bool is_valid (NT tol = NT(10) * std::numeric_limits<NT>::epsilon()) const; + + // POST: returns the time in seconds taken by the constructor call for + // computing the smallest enclosing ball + double get_time() const; + + // POST: deletes dynamically allocated arrays + ~Miniball(); + + private: + void mtf_mb (Sit n); + void mtf_move_to_front (Sit j); + void pivot_mb (Pit n); + void pivot_move_to_front (Pit j); + NT excess (Pit pit) const; + void pop (); + bool push (Pit pit); + NT suboptimality () const; + void create_arrays(); + void delete_arrays(); + }; + + // Class Definition + // ================ + template <typename CoordAccessor> + Miniball<CoordAccessor>::Miniball (int d_, Pit begin, Pit end, + CoordAccessor ca) + : d (d_), + points_begin (begin), + points_end (end), + coord_accessor (ca), + time (clock()), + nt0 (NT(0)), + L(), + support_end (L.begin()), + fsize(0), + ssize(0), + current_c (NULL), + current_sqr_r (NT(-1)), + c (NULL), + sqr_r (NULL), + q0 (NULL), + z (NULL), + f (NULL), + v (NULL), + a (NULL) + { + assert (points_begin != points_end); + create_arrays(); + + // set initial center + for (int j=0; j<d; ++j) c[0][j] = nt0; + current_c = c[0]; + + // compute miniball + pivot_mb (points_end); + + // update time + time = (clock() - time) / CLOCKS_PER_SEC; + } + + template <typename CoordAccessor> + Miniball<CoordAccessor>::~Miniball() + { + delete_arrays(); + } + + template <typename CoordAccessor> + void Miniball<CoordAccessor>::create_arrays() + { + c = new NT*[d+1]; + v = new NT*[d+1]; + a = new NT*[d+1]; + for (int i=0; i<d+1; ++i) { + c[i] = new NT[d]; + v[i] = new NT[d]; + a[i] = new NT[d]; + } + sqr_r = new NT[d+1]; + q0 = new NT[d]; + z = new NT[d+1]; + f = new NT[d+1]; + } + + template <typename CoordAccessor> + void Miniball<CoordAccessor>::delete_arrays() + { + delete[] f; + delete[] z; + delete[] q0; + delete[] sqr_r; + for (int i=0; i<d+1; ++i) { + delete[] a[i]; + delete[] v[i]; + delete[] c[i]; + } + delete[] a; + delete[] v; + delete[] c; + } + + template <typename CoordAccessor> + const typename Miniball<CoordAccessor>::NT* + Miniball<CoordAccessor>::center () const + { + return current_c; + } + + template <typename CoordAccessor> + typename Miniball<CoordAccessor>::NT + Miniball<CoordAccessor>::squared_radius () const + { + return current_sqr_r; + } + + template <typename CoordAccessor> + int Miniball<CoordAccessor>::nr_support_points () const + { + assert (ssize < d+2); + return ssize; + } + + template <typename CoordAccessor> + typename Miniball<CoordAccessor>::SupportPointIterator + Miniball<CoordAccessor>::support_points_begin () const + { + return L.begin(); + } + + template <typename CoordAccessor> + typename Miniball<CoordAccessor>::SupportPointIterator + Miniball<CoordAccessor>::support_points_end () const + { + return support_end; + } + + template <typename CoordAccessor> + typename Miniball<CoordAccessor>::NT + Miniball<CoordAccessor>::relative_error (NT& subopt) const + { + NT e, max_e = nt0; + // compute maximum absolute excess of support points + for (SupportPointIterator it = support_points_begin(); + it != support_points_end(); ++it) { + e = excess (*it); + if (e < nt0) e = -e; + if (e > max_e) { + max_e = e; + } + } + // compute maximum excess of any point + for (Pit i = points_begin; i != points_end; ++i) + if ((e = excess (i)) > max_e) + max_e = e; + + subopt = suboptimality(); + assert (current_sqr_r > nt0 || max_e == nt0); + return (current_sqr_r == nt0 ? nt0 : max_e / current_sqr_r); + } + + template <typename CoordAccessor> + bool Miniball<CoordAccessor>::is_valid (NT tol) const + { + NT suboptimality; + return ( (relative_error (suboptimality) <= tol) && (suboptimality == 0) ); + } + + template <typename CoordAccessor> + double Miniball<CoordAccessor>::get_time() const + { + return time; + } + + template <typename CoordAccessor> + void Miniball<CoordAccessor>::mtf_mb (Sit n) + { + // Algorithm 1: mtf_mb (L_{n-1}, B), where L_{n-1} = [L.begin, n) + // B: the set of forced points, defining the current ball + // S: the superset of support points computed by the algorithm + // -------------------------------------------------------------- + // from B. Gaertner, Fast and Robust Smallest Enclosing Balls, ESA 1999, + // http://www.inf.ethz.ch/personal/gaertner/texts/own_work/esa99_final.pdf + + // PRE: B = S + assert (fsize == ssize); + + support_end = L.begin(); + if ((fsize) == d+1) return; + + // incremental construction + for (Sit i = L.begin(); i != n;) + { + // INV: (support_end - L.begin() == |S|-|B|) + assert (std::distance (L.begin(), support_end) == ssize - fsize); + + Sit j = i++; + if (excess(*j) > nt0) + if (push(*j)) { // B := B + p_i + mtf_mb (j); // mtf_mb (L_{i-1}, B + p_i) + pop(); // B := B - p_i + mtf_move_to_front(j); + } + } + // POST: the range [L.begin(), support_end) stores the set S\B + } + + template <typename CoordAccessor> + void Miniball<CoordAccessor>::mtf_move_to_front (Sit j) + { + if (support_end == j) + support_end++; + L.splice (L.begin(), L, j); + } + + template <typename CoordAccessor> + void Miniball<CoordAccessor>::pivot_mb (Pit n) + { + // Algorithm 2: pivot_mb (L_{n-1}), where L_{n-1} = [L.begin, n) + // -------------------------------------------------------------- + // from B. Gaertner, Fast and Robust Smallest Enclosing Balls, ESA 1999, + // http://www.inf.ethz.ch/personal/gaertner/texts/own_work/esa99_final.pdf + NT old_sqr_r; + const NT* c; + Pit pivot, k; + NT e, max_e, sqr_r; + Cit p; + do { + old_sqr_r = current_sqr_r; + sqr_r = current_sqr_r; + + pivot = points_begin; + max_e = nt0; + for (k = points_begin; k != n; ++k) { + p = coord_accessor(k); + e = -sqr_r; + c = current_c; + for (int j=0; j<d; ++j) + e += mb_sqr<NT>(*p++-*c++); + if (e > max_e) { + max_e = e; + pivot = k; + } + } + + if (max_e > nt0) { + // check if the pivot is already contained in the support set + if (std::find(L.begin(), support_end, pivot) == support_end) { + assert (fsize == 0); + if (push (pivot)) { + mtf_mb(support_end); + pop(); + pivot_move_to_front(pivot); + } + } + } + } while (old_sqr_r < current_sqr_r); + } + + template <typename CoordAccessor> + void Miniball<CoordAccessor>::pivot_move_to_front (Pit j) + { + L.push_front(j); + if (std::distance(L.begin(), support_end) == d+2) + support_end--; + } + + template <typename CoordAccessor> + inline typename Miniball<CoordAccessor>::NT + Miniball<CoordAccessor>::excess (Pit pit) const + { + Cit p = coord_accessor(pit); + NT e = -current_sqr_r; + NT* c = current_c; + for (int k=0; k<d; ++k){ + e += mb_sqr<NT>(*p++-*c++); + } + return e; + } + + template <typename CoordAccessor> + void Miniball<CoordAccessor>::pop () + { + --fsize; + } + + template <typename CoordAccessor> + bool Miniball<CoordAccessor>::push (Pit pit) + { + int i, j; + NT eps = mb_sqr<NT>(std::numeric_limits<NT>::epsilon()); + + Cit cit = coord_accessor(pit); + Cit p = cit; + + if (fsize==0) { + for (i=0; i<d; ++i) + q0[i] = *p++; + for (i=0; i<d; ++i) + c[0][i] = q0[i]; + sqr_r[0] = nt0; + } + else { + // set v_fsize to Q_fsize + for (i=0; i<d; ++i) + //v[fsize][i] = p[i]-q0[i]; + v[fsize][i] = *p++-q0[i]; + + // compute the a_{fsize,i}, i< fsize + for (i=1; i<fsize; ++i) { + a[fsize][i] = nt0; + for (j=0; j<d; ++j) + a[fsize][i] += v[i][j] * v[fsize][j]; + a[fsize][i]*=(2/z[i]); + } + + // update v_fsize to Q_fsize-\bar{Q}_fsize + for (i=1; i<fsize; ++i) { + for (j=0; j<d; ++j) + v[fsize][j] -= a[fsize][i]*v[i][j]; + } + + // compute z_fsize + z[fsize]=nt0; + for (j=0; j<d; ++j) + z[fsize] += mb_sqr<NT>(v[fsize][j]); + z[fsize]*=2; + + // reject push if z_fsize too small + if (z[fsize]<eps*current_sqr_r) { + return false; + } + + // update c, sqr_r + p=cit; + NT e = -sqr_r[fsize-1]; + for (i=0; i<d; ++i) + e += mb_sqr<NT>(*p++-c[fsize-1][i]); + f[fsize]=e/z[fsize]; + + for (i=0; i<d; ++i) + c[fsize][i] = c[fsize-1][i]+f[fsize]*v[fsize][i]; + sqr_r[fsize] = sqr_r[fsize-1] + e*f[fsize]/2; + } + current_c = c[fsize]; + current_sqr_r = sqr_r[fsize]; + ssize = ++fsize; + return true; + } + + template <typename CoordAccessor> + typename Miniball<CoordAccessor>::NT + Miniball<CoordAccessor>::suboptimality () const + { + NT* l = new NT[d+1]; + NT min_l = nt0; + l[0] = NT(1); + for (int i=ssize-1; i>0; --i) { + l[i] = f[i]; + for (int k=ssize-1; k>i; --k) + l[i]-=a[k][i]*l[k]; + if (l[i] < min_l) min_l = l[i]; + l[0] -= l[i]; + } + if (l[0] < min_l) min_l = l[0]; + delete[] l; + if (min_l < nt0) + return -min_l; + return nt0; + } +} // namespace Miniball + +} // namespace Gudhi + +#endif // MINIBALL_HPP_ diff --git a/include/gudhi/Neighbors_finder.h b/include/gudhi/Neighbors_finder.h index 87c7cee5..36a63ea0 100644 --- a/include/gudhi/Neighbors_finder.h +++ b/include/gudhi/Neighbors_finder.h @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA + * Copyright (C) 2015 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 diff --git a/include/gudhi/Null_output_iterator.h b/include/gudhi/Null_output_iterator.h index 42e6e449..c700af5f 100644 --- a/include/gudhi/Null_output_iterator.h +++ b/include/gudhi/Null_output_iterator.h @@ -4,7 +4,7 @@ * * Author(s): Marc Glisse * - * Copyright (C) 2017 INRIA + * Copyright (C) 2017 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 diff --git a/include/gudhi/Off_reader.h b/include/gudhi/Off_reader.h index 4fcd2af2..05a1e145 100644 --- a/include/gudhi/Off_reader.h +++ b/include/gudhi/Off_reader.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 @@ -105,25 +105,26 @@ class Off_reader { bool is_off_file = (line.find("OFF") != std::string::npos); bool is_noff_file = (line.find("nOFF") != std::string::npos); + + if (!is_off_file && !is_noff_file) { std::cerr << line << std::endl; std::cerr << "missing off header\n"; return false; } + if (is_noff_file) { + // Should be on a separate line, but we accept it on the same line as the number of vertices + stream_ >> off_info_.dim; + } else { + off_info_.dim = 3; + } + if (!goto_next_uncomment_line(line)) return false; std::istringstream iss(line); - if ((is_off_file) && (!is_noff_file)) { - off_info_.dim = 3; - if (!(iss >> off_info_.num_vertices >> off_info_.num_faces >> off_info_.num_edges)) { - std::cerr << "incorrect number of vertices/faces/edges\n"; - return false; - } - } else { - if (!(iss >> off_info_.dim >> off_info_.num_vertices >> off_info_.num_faces >> off_info_.num_edges)) { + if (!(iss >> off_info_.num_vertices >> off_info_.num_faces >> off_info_.num_edges)) { std::cerr << "incorrect number of vertices/faces/edges\n"; return false; - } } off_visitor.init(off_info_.dim, off_info_.num_vertices, off_info_.num_faces, off_info_.num_edges); @@ -131,10 +132,12 @@ class Off_reader { } bool goto_next_uncomment_line(std::string& uncomment_line) { - uncomment_line.clear(); - do - std::getline(stream_, uncomment_line); while (uncomment_line[0] == '%'); - return (uncomment_line.size() > 0 && uncomment_line[0] != '%'); + do { + // skip whitespace, including empty lines + if (!std::ifstream::sentry(stream_)) return false; + std::getline(stream_, uncomment_line); + } while (uncomment_line[0] == '#'); + return static_cast<bool>(stream_); } template<typename OffVisitor> diff --git a/include/gudhi/PSSK.h b/include/gudhi/PSSK.h index 630f5623..e1174455 100644 --- a/include/gudhi/PSSK.h +++ b/include/gudhi/PSSK.h @@ -4,7 +4,7 @@ * * Author(s): Pawel Dlotko * - * Copyright (C) 2016 INRIA (France) + * 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 diff --git a/include/gudhi/Persistence_graph.h b/include/gudhi/Persistence_graph.h index 622b0691..cb163623 100644 --- a/include/gudhi/Persistence_graph.h +++ b/include/gudhi/Persistence_graph.h @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA + * Copyright (C) 2015 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 diff --git a/include/gudhi/Persistence_heat_maps.h b/include/gudhi/Persistence_heat_maps.h index a80c3c40..35e51e63 100644 --- a/include/gudhi/Persistence_heat_maps.h +++ b/include/gudhi/Persistence_heat_maps.h @@ -4,7 +4,7 @@ * * Author(s): Pawel Dlotko * - * Copyright (C) 2016 INRIA (France) + * 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 diff --git a/include/gudhi/Persistence_intervals.h b/include/gudhi/Persistence_intervals.h index 3d04d8b7..76eac7d7 100644 --- a/include/gudhi/Persistence_intervals.h +++ b/include/gudhi/Persistence_intervals.h @@ -4,7 +4,7 @@ * * Author(s): Pawel Dlotko * - * Copyright (C) 2016 INRIA (France) + * 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 diff --git a/include/gudhi/Persistence_intervals_with_distances.h b/include/gudhi/Persistence_intervals_with_distances.h index 79908883..f48d1a3b 100644 --- a/include/gudhi/Persistence_intervals_with_distances.h +++ b/include/gudhi/Persistence_intervals_with_distances.h @@ -4,7 +4,7 @@ * * Author(s): Pawel Dlotko * - * Copyright (C) 2016 INRIA (France) + * 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 diff --git a/include/gudhi/Persistence_landscape.h b/include/gudhi/Persistence_landscape.h index c5aa7867..9cab0166 100644 --- a/include/gudhi/Persistence_landscape.h +++ b/include/gudhi/Persistence_landscape.h @@ -4,7 +4,7 @@ * * Author(s): Pawel Dlotko * - * Copyright (C) 2016 INRIA (France) + * 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 @@ -734,7 +734,7 @@ double Persistence_landscape::compute_integral_of_landscape(double p) const { double Persistence_landscape::compute_value_at_a_given_point(unsigned level, double x) const { bool compute_value_at_a_given_pointDbg = false; // in such a case lambda_level = 0. - if (level > this->land.size()) return 0; + if (level >= this->land.size()) return 0; // we know that the points in this->land[level] are ordered according to x coordinate. Therefore, we can find the // point by using bisection: @@ -1235,40 +1235,43 @@ double compute_inner_product(const Persistence_landscape& l1, const Persistence_ std::cerr << "Computing inner product for a level : " << level << std::endl; getchar(); } - if (l1.land[level].size() * l2.land[level].size() == 0) continue; + auto&& l1_land_level = l1.land[level]; + auto&& l2_land_level = l2.land[level]; + + if (l1_land_level.size() * l2_land_level.size() == 0) continue; // endpoints of the interval on which we will compute the inner product of two locally linear functions: double x1 = -std::numeric_limits<int>::max(); double x2; - if (l1.land[level][1].first < l2.land[level][1].first) { - x2 = l1.land[level][1].first; + if (l1_land_level[1].first < l2_land_level[1].first) { + x2 = l1_land_level[1].first; } else { - x2 = l2.land[level][1].first; + x2 = l2_land_level[1].first; } // iterators for the landscapes l1 and l2 size_t l1It = 0; size_t l2It = 0; - while ((l1It < l1.land[level].size() - 1) && (l2It < l2.land[level].size() - 1)) { + while ((l1It < l1_land_level.size() - 1) && (l2It < l2_land_level.size() - 1)) { // compute the value of a inner product on a interval [x1,x2] double a, b, c, d; - if (l1.land[level][l1It + 1].first != l1.land[level][l1It].first) { - a = (l1.land[level][l1It + 1].second - l1.land[level][l1It].second) / - (l1.land[level][l1It + 1].first - l1.land[level][l1It].first); + if (l1_land_level[l1It + 1].first != l1_land_level[l1It].first) { + a = (l1_land_level[l1It + 1].second - l1_land_level[l1It].second) / + (l1_land_level[l1It + 1].first - l1_land_level[l1It].first); } else { a = 0; } - b = l1.land[level][l1It].second - a * l1.land[level][l1It].first; - if (l2.land[level][l2It + 1].first != l2.land[level][l2It].first) { - c = (l2.land[level][l2It + 1].second - l2.land[level][l2It].second) / - (l2.land[level][l2It + 1].first - l2.land[level][l2It].first); + b = l1_land_level[l1It].second - a * l1_land_level[l1It].first; + if (l2_land_level[l2It + 1].first != l2_land_level[l2It].first) { + c = (l2_land_level[l2It + 1].second - l2_land_level[l2It].second) / + (l2_land_level[l2It + 1].first - l2_land_level[l2It].first); } else { c = 0; } - d = l2.land[level][l2It].second - c * l2.land[level][l2It].first; + d = l2_land_level[l2It].second - c * l2_land_level[l2It].first; double contributionFromThisPart = (a * c * x2 * x2 * x2 / 3 + (a * d + b * c) * x2 * x2 / 2 + b * d * x2) - (a * c * x1 * x1 * x1 / 3 + (a * d + b * c) * x1 * x1 / 2 + b * d * x1); @@ -1276,10 +1279,10 @@ double compute_inner_product(const Persistence_landscape& l1, const Persistence_ result += contributionFromThisPart; if (dbg) { - std::cerr << "[l1.land[level][l1It].first,l1.land[level][l1It+1].first] : " << l1.land[level][l1It].first - << " , " << l1.land[level][l1It + 1].first << std::endl; - std::cerr << "[l2.land[level][l2It].first,l2.land[level][l2It+1].first] : " << l2.land[level][l2It].first - << " , " << l2.land[level][l2It + 1].first << std::endl; + std::cerr << "[l1_land_level[l1It].first,l1_land_level[l1It+1].first] : " << l1_land_level[l1It].first + << " , " << l1_land_level[l1It + 1].first << std::endl; + std::cerr << "[l2_land_level[l2It].first,l2_land_level[l2It+1].first] : " << l2_land_level[l2It].first + << " , " << l2_land_level[l2It + 1].first << std::endl; std::cerr << "a : " << a << ", b : " << b << " , c: " << c << ", d : " << d << std::endl; std::cerr << "x1 : " << x1 << " , x2 : " << x2 << std::endl; std::cerr << "contributionFromThisPart : " << contributionFromThisPart << std::endl; @@ -1288,14 +1291,14 @@ double compute_inner_product(const Persistence_landscape& l1, const Persistence_ } // we have two intervals in which functions are constant: - // [l1.land[level][l1It].first , l1.land[level][l1It+1].first] + // [l1_land_level[l1It].first , l1_land_level[l1It+1].first] // and - // [l2.land[level][l2It].first , l2.land[level][l2It+1].first] + // [l2_land_level[l2It].first , l2_land_level[l2It+1].first] // We also have an interval [x1,x2]. Since the intervals in the landscapes cover the whole R, then it is clear // that x2 - // is either l1.land[level][l1It+1].first of l2.land[level][l2It+1].first or both. Lets test it. - if (x2 == l1.land[level][l1It + 1].first) { - if (x2 == l2.land[level][l2It + 1].first) { + // is either l1_land_level[l1It+1].first of l2_land_level[l2It+1].first or both. Lets test it. + if (x2 == l1_land_level[l1It + 1].first) { + if (x2 == l2_land_level[l2It + 1].first) { // in this case, we increment both: ++l2It; if (dbg) { @@ -1314,12 +1317,16 @@ double compute_inner_product(const Persistence_landscape& l1, const Persistence_ std::cerr << "Incrementing second \n"; } } + + if ( l1It + 1 >= l1_land_level.size() )break; + if ( l2It + 1 >= l2_land_level.size() )break; + // Now, we shift x1 and x2: x1 = x2; - if (l1.land[level][l1It + 1].first < l2.land[level][l2It + 1].first) { - x2 = l1.land[level][l1It + 1].first; + if (l1_land_level[l1It + 1].first < l2_land_level[l2It + 1].first) { + x2 = l1_land_level[l1It + 1].first; } else { - x2 = l2.land[level][l2It + 1].first; + x2 = l2_land_level[l2It + 1].first; } } } diff --git a/include/gudhi/Persistence_landscape_on_grid.h b/include/gudhi/Persistence_landscape_on_grid.h index 84fd22ed..fd8a181c 100644 --- a/include/gudhi/Persistence_landscape_on_grid.h +++ b/include/gudhi/Persistence_landscape_on_grid.h @@ -4,7 +4,7 @@ * * Author(s): Pawel Dlotko * - * Copyright (C) 2016 INRIA (France) + * 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 diff --git a/include/gudhi/Persistence_vectors.h b/include/gudhi/Persistence_vectors.h index 63577e46..9c04be1d 100644 --- a/include/gudhi/Persistence_vectors.h +++ b/include/gudhi/Persistence_vectors.h @@ -4,7 +4,7 @@ * * Author(s): Pawel Dlotko * - * Copyright (C) 2016 INRIA (France) + * 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 diff --git a/include/gudhi/Persistent_cohomology.h b/include/gudhi/Persistent_cohomology.h index e0a147b3..c68b5c0b 100644 --- a/include/gudhi/Persistent_cohomology.h +++ b/include/gudhi/Persistent_cohomology.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 @@ -285,7 +285,7 @@ class Persistent_cohomology { } } cpx_->assign_key(sigma, cpx_->null_key()); - } else { // If ku == kv, same connected component: create a 1-cocycle class. + } else if (dim_max_ > 1) { // If ku == kv, same connected component: create a 1-cocycle class. create_cocycle(sigma, coeff_field_.multiplicative_identity(), coeff_field_.characteristic()); } } diff --git a/include/gudhi/Persistent_cohomology/Field_Zp.h b/include/gudhi/Persistent_cohomology/Field_Zp.h index 6db16e69..e98b4bb4 100644 --- a/include/gudhi/Persistent_cohomology/Field_Zp.h +++ b/include/gudhi/Persistent_cohomology/Field_Zp.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 diff --git a/include/gudhi/Persistent_cohomology/Multi_field.h b/include/gudhi/Persistent_cohomology/Multi_field.h index 38bc08d1..2bae8654 100644 --- a/include/gudhi/Persistent_cohomology/Multi_field.h +++ b/include/gudhi/Persistent_cohomology/Multi_field.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 diff --git a/include/gudhi/Persistent_cohomology/Persistent_cohomology_column.h b/include/gudhi/Persistent_cohomology/Persistent_cohomology_column.h index 5deb2d88..de6c0750 100644 --- a/include/gudhi/Persistent_cohomology/Persistent_cohomology_column.h +++ b/include/gudhi/Persistent_cohomology/Persistent_cohomology_column.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 diff --git a/include/gudhi/Point.h b/include/gudhi/Point.h index 0479e71e..345a8465 100644 --- a/include/gudhi/Point.h +++ b/include/gudhi/Point.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 diff --git a/include/gudhi/Points_3D_off_io.h b/include/gudhi/Points_3D_off_io.h index b0d24998..704f73a7 100644 --- a/include/gudhi/Points_3D_off_io.h +++ b/include/gudhi/Points_3D_off_io.h @@ -4,7 +4,7 @@ * * Author(s): Vincent Rouvreau * - * Copyright (C) 2015 INRIA Saclay (France) + * Copyright (C) 2015 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 diff --git a/include/gudhi/Points_off_io.h b/include/gudhi/Points_off_io.h index 29af8a8a..38029658 100644 --- a/include/gudhi/Points_off_io.h +++ b/include/gudhi/Points_off_io.h @@ -4,7 +4,7 @@ * * Author(s): Vincent Rouvreau * - * Copyright (C) 2015 INRIA Saclay (France) + * Copyright (C) 2015 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 @@ -126,9 +126,9 @@ class Points_off_visitor_reader { * \code $> ./vector_double_off_reader ../../data/points/alphacomplexdoc.off * \endcode * - * the program output is: + * the program outputs a file ../../data/points/alphacomplexdoc.off.txt: * - * \include common/cgaloffreader_result.txt + * \include common/vectordoubleoffreader_result.txt */ template<typename Point_d> class Points_off_reader { diff --git a/include/gudhi/Rips_complex.h b/include/gudhi/Rips_complex.h index 1e4b76a7..f0fe57f4 100644 --- a/include/gudhi/Rips_complex.h +++ b/include/gudhi/Rips_complex.h @@ -4,7 +4,7 @@ * * Author(s): Clément Maria, Pawel Dlotko, Vincent Rouvreau * - * Copyright (C) 2016 INRIA + * 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 diff --git a/include/gudhi/Simple_object_pool.h b/include/gudhi/Simple_object_pool.h index fb9c8e23..47283521 100644 --- a/include/gudhi/Simple_object_pool.h +++ b/include/gudhi/Simple_object_pool.h @@ -4,7 +4,7 @@ * * Author(s): Marc Glisse * - * Copyright (C) 2015 INRIA Saclay - Ile de France + * Copyright (C) 2015 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 diff --git a/include/gudhi/Simplex_tree.h b/include/gudhi/Simplex_tree.h index 7456cb1f..ee96d5a2 100644 --- a/include/gudhi/Simplex_tree.h +++ b/include/gudhi/Simplex_tree.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 @@ -689,7 +689,11 @@ class Simplex_tree { return { null_simplex(), true }; // ----->> // Copy before sorting - thread_local std::vector<Vertex_handle> copy; + // Thread local is not available on XCode version < V.8 - It will slow down computation +#ifdef GUDHI_CAN_USE_CXX11_THREAD_LOCAL + thread_local +#endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL + std::vector<Vertex_handle> copy; copy.clear(); copy.insert(copy.end(), first, last); std::sort(std::begin(copy), std::end(copy)); @@ -1238,9 +1242,8 @@ class Simplex_tree { } public: - /** \brief Browse the simplex tree to ensure the filtration is not decreasing. - * The simplex tree is browsed starting from the root until the leaf, and the filtration values are set with their - * parent value (increased), in case the values are decreasing. + /** \brief This function ensures that each simplex has a higher filtration value than its faces by increasing the + * filtration values. * @return The filtration modification information. * \post Some simplex tree functions require the filtration to be valid. `make_filtration_non_decreasing()` * function is not launching `initialize_filtration()` but returns the filtration modification information. If the diff --git a/include/gudhi/Simplex_tree/Simplex_tree_iterators.h b/include/gudhi/Simplex_tree/Simplex_tree_iterators.h index ab7346d4..02c8bb64 100644 --- a/include/gudhi/Simplex_tree/Simplex_tree_iterators.h +++ b/include/gudhi/Simplex_tree/Simplex_tree_iterators.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 @@ -101,7 +101,9 @@ class Simplex_tree_boundary_simplex_iterator : public boost::iterator_facade< // any end() iterator explicit Simplex_tree_boundary_simplex_iterator(SimplexTree * st) - : sib_(nullptr), + : last_(st->null_vertex()), + next_(st->null_vertex()), + sib_(nullptr), sh_(st->null_simplex()), st_(st) { } @@ -109,7 +111,9 @@ class Simplex_tree_boundary_simplex_iterator : public boost::iterator_facade< template<class SimplexHandle> Simplex_tree_boundary_simplex_iterator(SimplexTree * st, SimplexHandle sh) : last_(sh->first), + next_(st->null_vertex()), sib_(nullptr), + sh_(st->null_simplex()), st_(st) { // Only check once at the beginning instead of for every increment, as this is expensive. if (SimplexTree::Options::contiguous_vertices) @@ -123,9 +127,7 @@ class Simplex_tree_boundary_simplex_iterator : public boost::iterator_facade< sh_ = sib_->members_.begin()+next_; else sh_ = sib_->find(next_); - } else { - sh_ = st->null_simplex(); - } // vertex: == end() + } } private: diff --git a/include/gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h b/include/gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h index 25d4888a..3a75ec72 100644 --- a/include/gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h +++ b/include/gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.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 diff --git a/include/gudhi/Simplex_tree/Simplex_tree_siblings.h b/include/gudhi/Simplex_tree/Simplex_tree_siblings.h index 1eca7f6f..ab2ca707 100644 --- a/include/gudhi/Simplex_tree/Simplex_tree_siblings.h +++ b/include/gudhi/Simplex_tree/Simplex_tree_siblings.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 diff --git a/include/gudhi/Simplex_tree/indexing_tag.h b/include/gudhi/Simplex_tree/indexing_tag.h index 0adeb46d..ec4461f3 100644 --- a/include/gudhi/Simplex_tree/indexing_tag.h +++ b/include/gudhi/Simplex_tree/indexing_tag.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 diff --git a/include/gudhi/Skeleton_blocker.h b/include/gudhi/Skeleton_blocker.h index aca2aa57..e8b6fde8 100644 --- a/include/gudhi/Skeleton_blocker.h +++ b/include/gudhi/Skeleton_blocker.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h b/include/gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h index ba3636bc..6c6a8638 100644 --- a/include/gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h +++ b/include/gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h b/include/gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h index d4b60613..feab7b3f 100644 --- a/include/gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h +++ b/include/gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h b/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h index 747e60f1..56009daf 100644 --- a/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h +++ b/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_geometric_traits.h b/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_geometric_traits.h index 275376e6..22c1668e 100644 --- a/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_geometric_traits.h +++ b/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_geometric_traits.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h b/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h index 3835cf77..144f1fd0 100644 --- a/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h +++ b/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h b/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h index aa6f2215..d7193157 100644 --- a/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h +++ b/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h b/include/gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h index fadf6619..dbfb4042 100644 --- a/include/gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h +++ b/include/gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/Skeleton_blocker/internal/Top_faces.h b/include/gudhi/Skeleton_blocker/internal/Top_faces.h index 2b681752..f80ca4fe 100644 --- a/include/gudhi/Skeleton_blocker/internal/Top_faces.h +++ b/include/gudhi/Skeleton_blocker/internal/Top_faces.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/Skeleton_blocker/internal/Trie.h b/include/gudhi/Skeleton_blocker/internal/Trie.h index 2c9602fa..7a5d38eb 100644 --- a/include/gudhi/Skeleton_blocker/internal/Trie.h +++ b/include/gudhi/Skeleton_blocker/internal/Trie.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_blockers_iterators.h b/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_blockers_iterators.h index d2fff960..95c5f7ef 100644 --- a/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_blockers_iterators.h +++ b/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_blockers_iterators.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_edges_iterators.h b/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_edges_iterators.h index b90dcf34..5c725aae 100644 --- a/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_edges_iterators.h +++ b/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_edges_iterators.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_iterators.h b/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_iterators.h index 1351614f..8054e64f 100644 --- a/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_iterators.h +++ b/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_iterators.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h b/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h index 2acdb555..e2024652 100644 --- a/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h +++ b/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_triangles_iterators.h b/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_triangles_iterators.h index 736941dd..a834fe1d 100644 --- a/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_triangles_iterators.h +++ b/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_triangles_iterators.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_vertices_iterators.h b/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_vertices_iterators.h index 9e9ae961..3a638ae6 100644 --- a/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_vertices_iterators.h +++ b/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_vertices_iterators.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/Skeleton_blocker_complex.h b/include/gudhi/Skeleton_blocker_complex.h index 4f052ba5..addd8104 100644 --- a/include/gudhi/Skeleton_blocker_complex.h +++ b/include/gudhi/Skeleton_blocker_complex.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/Skeleton_blocker_contractor.h b/include/gudhi/Skeleton_blocker_contractor.h index df884c93..13086161 100644 --- a/include/gudhi/Skeleton_blocker_contractor.h +++ b/include/gudhi/Skeleton_blocker_contractor.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA Sophia Antipolis-Mediterranee (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 diff --git a/include/gudhi/Skeleton_blocker_geometric_complex.h b/include/gudhi/Skeleton_blocker_geometric_complex.h index 95331b7a..39b88ceb 100644 --- a/include/gudhi/Skeleton_blocker_geometric_complex.h +++ b/include/gudhi/Skeleton_blocker_geometric_complex.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/Skeleton_blocker_link_complex.h b/include/gudhi/Skeleton_blocker_link_complex.h index 4db075b0..428d4e9b 100644 --- a/include/gudhi/Skeleton_blocker_link_complex.h +++ b/include/gudhi/Skeleton_blocker_link_complex.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/Skeleton_blocker_simplifiable_complex.h b/include/gudhi/Skeleton_blocker_simplifiable_complex.h index 544e02e8..d5adb39d 100644 --- a/include/gudhi/Skeleton_blocker_simplifiable_complex.h +++ b/include/gudhi/Skeleton_blocker_simplifiable_complex.h @@ -4,7 +4,7 @@ * * Author(s): David Salinas * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/Sparse_rips_complex.h b/include/gudhi/Sparse_rips_complex.h new file mode 100644 index 00000000..4dcc08ed --- /dev/null +++ b/include/gudhi/Sparse_rips_complex.h @@ -0,0 +1,175 @@ +/* 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): Marc Glisse + * + * Copyright (C) 2018 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 SPARSE_RIPS_COMPLEX_H_ +#define SPARSE_RIPS_COMPLEX_H_ + +#include <gudhi/Debug_utils.h> +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/choose_n_farthest_points.h> + +#include <boost/graph/adjacency_list.hpp> +#include <boost/range/metafunctions.hpp> + +#include <vector> + +namespace Gudhi { + +namespace rips_complex { + +// The whole interface is copied on Rips_complex. A redesign should be discussed with all complex creation classes in +// mind. + +/** + * \class Sparse_rips_complex + * \brief Sparse Rips complex data structure. + * + * \ingroup rips_complex + * + * \details + * This class is used to construct a sparse \f$(1+O(\epsilon))\f$-approximation of `Rips_complex`, i.e. a filtered + * simplicial complex that is multiplicatively \f$(1+O(\epsilon))\f$-interleaved with the Rips filtration. + * + * \tparam Filtration_value is the type used to store the filtration values of the simplicial complex. + */ +template <typename Filtration_value> +class Sparse_rips_complex { + private: + // TODO(MG): use a different graph where we know we can safely insert in parallel. + typedef typename boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS, + boost::property<vertex_filtration_t, Filtration_value>, + boost::property<edge_filtration_t, Filtration_value>> + Graph; + + typedef int Vertex_handle; + + public: + /** \brief Sparse_rips_complex constructor from a list of points. + * + * @param[in] points Range of points. + * @param[in] distance Distance function that returns a `Filtration_value` from 2 given points. + * @param[in] epsilon Approximation parameter. epsilon must be positive. + * + */ + template <typename RandomAccessPointRange, typename Distance> + Sparse_rips_complex(const RandomAccessPointRange& points, Distance distance, double epsilon) { + GUDHI_CHECK(epsilon > 0, "epsilon must be positive"); + std::vector<Vertex_handle> sorted_points; + std::vector<Filtration_value> params; + auto dist_fun = [&](Vertex_handle i, Vertex_handle j) { return distance(points[i], points[j]); }; + Ker<decltype(dist_fun)> kernel(dist_fun); + subsampling::choose_n_farthest_points(kernel, boost::irange<Vertex_handle>(0, boost::size(points)), -1, -1, + std::back_inserter(sorted_points), std::back_inserter(params)); + compute_sparse_graph(sorted_points, params, dist_fun, epsilon); + } + + /** \brief Sparse_rips_complex constructor from a distance matrix. + * + * @param[in] distance_matrix Range of range of distances. + * `distance_matrix[i][j]` returns the distance between points \f$i\f$ and + * \f$j\f$ as long as \f$ 0 \leqslant i < j \leqslant + * distance\_matrix.size().\f$ + * @param[in] epsilon Approximation parameter. epsilon must be positive. + */ + template <typename DistanceMatrix> + Sparse_rips_complex(const DistanceMatrix& distance_matrix, double epsilon) + : Sparse_rips_complex(boost::irange<Vertex_handle>(0, boost::size(distance_matrix)), + [&](Vertex_handle i, Vertex_handle j) { return distance_matrix[j][i]; }, epsilon) {} + + /** \brief Fills the simplicial complex with the sparse Rips graph and + * expands it with all the cliques, stopping at a given maximal dimension. + * + * \tparam SimplicialComplexForRips must meet `SimplicialComplexForRips` concept. + * + * @param[in] complex the complex to fill + * @param[in] dim_max maximal dimension of the simplicial complex. + * @exception std::invalid_argument In debug mode, if `complex.num_vertices()` does not return 0. + * + */ + template <typename SimplicialComplexForRips> + void create_complex(SimplicialComplexForRips& complex, int dim_max) { + GUDHI_CHECK(complex.num_vertices() == 0, + std::invalid_argument("Sparse_rips_complex::create_complex - simplicial complex is not empty")); + + complex.insert_graph(graph_); + complex.expansion(dim_max); + } + + private: + // choose_n_farthest_points wants the distance function in this form... + template <class Distance> + struct Ker { + typedef std::size_t Point_d; // index into point range + Ker(Distance& d) : dist(d) {} + // Despite the name, this is not squared... + typedef Distance Squared_distance_d; + Squared_distance_d& squared_distance_d_object() const { return dist; } + Distance& dist; + }; + + // PointRange must be random access. + template <typename PointRange, typename ParamRange, typename Distance> + void compute_sparse_graph(const PointRange& points, const ParamRange& params, Distance& dist, double epsilon) { + const int n = boost::size(points); + graph_.~Graph(); + new (&graph_) Graph(n); + // for(auto v : vertices(g)) // doesn't work :-( + typename boost::graph_traits<Graph>::vertex_iterator v_i, v_e; + for (std::tie(v_i, v_e) = vertices(graph_); v_i != v_e; ++v_i) { + auto v = *v_i; + // This whole loop might not be necessary, leave it until someone investigates if it is safe to remove. + put(vertex_filtration_t(), graph_, v, 0); + } + + // TODO(MG): + // - make it parallel + // - only test near-enough neighbors + for (int i = 0; i < n; ++i) + for (int j = i + 1; j < n; ++j) { + auto&& pi = points[i]; + auto&& pj = points[j]; + auto d = dist(pi, pj); + auto li = params[i]; + auto lj = params[j]; + GUDHI_CHECK(lj <= li, "Bad furthest point sorting"); + Filtration_value alpha; + + // The paper has d/2 and d-lj/e to match the Cech, but we use doubles to match the Rips + if (d * epsilon <= 2 * lj) + alpha = d; + else if (d * epsilon <= li + lj && (epsilon >= 1 || d * epsilon <= lj * (1 + 1 / (1 - epsilon)))) + alpha = (d - lj / epsilon) * 2; + else + continue; + + add_edge(pi, pj, alpha, graph_); + } + } + + Graph graph_; +}; + +} // namespace rips_complex + +} // namespace Gudhi + +#endif // SPARSE_RIPS_COMPLEX_H_ diff --git a/include/gudhi/Strong_witness_complex.h b/include/gudhi/Strong_witness_complex.h index b3d00b11..fd6b3f38 100644 --- a/include/gudhi/Strong_witness_complex.h +++ b/include/gudhi/Strong_witness_complex.h @@ -4,7 +4,7 @@ * * Author(s): Siargey Kachanovich * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 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 diff --git a/include/gudhi/Tangential_complex.h b/include/gudhi/Tangential_complex.h index 6f061922..9d8fdcd3 100644 --- a/include/gudhi/Tangential_complex.h +++ b/include/gudhi/Tangential_complex.h @@ -4,7 +4,7 @@ * * Author(s): Clement Jamin * - * Copyright (C) 2016 INRIA + * 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 @@ -1100,7 +1100,10 @@ class Tangential_complex { // of the sphere "star sphere" centered at "center_vertex" // and which contains all the // circumspheres of the star of "center_vertex" - boost::optional<FT> squared_star_sphere_radius_plus_margin; + boost::optional<FT> squared_star_sphere_radius_plus_margin = boost::make_optional(false, FT()); + // This is the strange way boost is recommending to get rid of "may be used uninitialized in this function". + // Former code was : + // boost::optional<FT> squared_star_sphere_radius_plus_margin; // Insert points until we find a point which is outside "star sphere" for (auto nn_it = ins_range.begin(); diff --git a/include/gudhi/Tangential_complex/Simplicial_complex.h b/include/gudhi/Tangential_complex/Simplicial_complex.h index 65c74ca5..f79186b0 100644 --- a/include/gudhi/Tangential_complex/Simplicial_complex.h +++ b/include/gudhi/Tangential_complex/Simplicial_complex.h @@ -4,7 +4,7 @@ * * Author(s): Clement Jamin * - * Copyright (C) 2016 INRIA + * 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 diff --git a/include/gudhi/Tangential_complex/config.h b/include/gudhi/Tangential_complex/config.h index ffefcd6b..e1af1ea6 100644 --- a/include/gudhi/Tangential_complex/config.h +++ b/include/gudhi/Tangential_complex/config.h @@ -4,7 +4,7 @@ * * Author(s): Clement Jamin * - * Copyright (C) 2016 INRIA + * 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 diff --git a/include/gudhi/Tangential_complex/utilities.h b/include/gudhi/Tangential_complex/utilities.h index b2d6d674..2dd46118 100644 --- a/include/gudhi/Tangential_complex/utilities.h +++ b/include/gudhi/Tangential_complex/utilities.h @@ -4,7 +4,7 @@ * * Author(s): Clement Jamin * - * Copyright (C) 2016 INRIA + * 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 diff --git a/include/gudhi/Unitary_tests_utils.h b/include/gudhi/Unitary_tests_utils.h index 8394a062..e07c8d42 100644 --- a/include/gudhi/Unitary_tests_utils.h +++ b/include/gudhi/Unitary_tests_utils.h @@ -4,7 +4,7 @@ * * Author(s): Vincent Rouvreau * - * Copyright (C) 2017 INRIA + * Copyright (C) 2017 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 diff --git a/include/gudhi/Witness_complex.h b/include/gudhi/Witness_complex.h index 53c38520..67885258 100644 --- a/include/gudhi/Witness_complex.h +++ b/include/gudhi/Witness_complex.h @@ -4,7 +4,7 @@ * * Author(s): Siargey Kachanovich * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 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 diff --git a/include/gudhi/Witness_complex/all_faces_in.h b/include/gudhi/Witness_complex/all_faces_in.h index b68d75a1..c7b732b9 100644 --- a/include/gudhi/Witness_complex/all_faces_in.h +++ b/include/gudhi/Witness_complex/all_faces_in.h @@ -4,7 +4,7 @@ * * Author(s): Siargey Kachanovich * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 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 diff --git a/include/gudhi/allocator.h b/include/gudhi/allocator.h index 4ede14e4..3de16a49 100644 --- a/include/gudhi/allocator.h +++ b/include/gudhi/allocator.h @@ -4,7 +4,7 @@ * * Author(s): Marc Glisse * - * Copyright (C) 2015 INRIA Saclay - Ile de France + * Copyright (C) 2015 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 diff --git a/include/gudhi/choose_n_farthest_points.h b/include/gudhi/choose_n_farthest_points.h index 8390b4c9..ab1c4c73 100644 --- a/include/gudhi/choose_n_farthest_points.h +++ b/include/gudhi/choose_n_farthest_points.h @@ -4,7 +4,7 @@ * * Author(s): Siargey Kachanovich * - * Copyright (C) 2016 INRIA + * 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 diff --git a/include/gudhi/common_persistence_representations.h b/include/gudhi/common_persistence_representations.h index 44e125a7..3d03f1f6 100644 --- a/include/gudhi/common_persistence_representations.h +++ b/include/gudhi/common_persistence_representations.h @@ -4,7 +4,7 @@ * * Author(s): Pawel Dlotko * - * Copyright (C) 2016 INRIA (France) + * 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 diff --git a/include/gudhi/console_color.h b/include/gudhi/console_color.h index c4671da3..a493e0d0 100644 --- a/include/gudhi/console_color.h +++ b/include/gudhi/console_color.h @@ -4,7 +4,7 @@ * * Author(s): Clement Jamin * - * Copyright (C) 2016 INRIA Sophia-Antipolis (France) + * 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 diff --git a/include/gudhi/distance_functions.h b/include/gudhi/distance_functions.h index 3a5d1fd5..5ef12f2e 100644 --- a/include/gudhi/distance_functions.h +++ b/include/gudhi/distance_functions.h @@ -4,7 +4,7 @@ * * Author(s): Clément Maria * - * Copyright (C) 2014 INRIA + * 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 @@ -25,7 +25,10 @@ #include <gudhi/Debug_utils.h> +#include <gudhi/Miniball.hpp> + #include <boost/range/metafunctions.hpp> +#include <boost/range/size.hpp> #include <cmath> // for std::sqrt #include <type_traits> // for std::decay @@ -68,6 +71,53 @@ class Euclidean_distance { } }; +/** @brief Compute the radius of the minimal enclosing ball between Points given by a range of coordinates. + * The points are assumed to have the same dimension. */ +class Minimal_enclosing_ball_radius { + public: + /** \brief Minimal_enclosing_ball_radius from two points. + * + * @param[in] point_1 First point. + * @param[in] point_2 second point. + * @return The minimal enclosing ball radius for the two points (aka. Euclidean distance / 2.). + * + * \tparam Point must be a range of Cartesian coordinates. + * + */ + template< typename Point > + typename std::iterator_traits<typename boost::range_iterator<Point>::type>::value_type + operator()(const Point& point_1, const Point& point_2) const { + return Euclidean_distance()(point_1, point_2) / 2.; + } + /** \brief Minimal_enclosing_ball_radius from a point cloud. + * + * @param[in] point_cloud The points. + * @return The minimal enclosing ball radius for the points. + * + * \tparam Point_cloud must be a range of points with Cartesian coordinates. + * Point_cloud is a range over a range of Coordinate. + * + */ + template< typename Point_cloud, + typename Point_iterator = typename boost::range_const_iterator<Point_cloud>::type, + typename Point = typename std::iterator_traits<Point_iterator>::value_type, + typename Coordinate_iterator = typename boost::range_const_iterator<Point>::type, + typename Coordinate = typename std::iterator_traits<Coordinate_iterator>::value_type> + Coordinate + operator()(const Point_cloud& point_cloud) const { + using Min_sphere = Miniball::Miniball<Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>; + + Min_sphere ms(boost::size(*point_cloud.begin()), point_cloud.begin(), point_cloud.end()); +#ifdef DEBUG_TRACES + std::cout << "Minimal_enclosing_ball_radius = " << std::sqrt(ms.squared_radius()) << " | nb points = " + << boost::size(point_cloud) << " | dimension = " + << boost::size(*point_cloud.begin()) << std::endl; +#endif // DEBUG_TRACES + + return std::sqrt(ms.squared_radius()); + } +}; + } // namespace Gudhi #endif // DISTANCE_FUNCTIONS_H_ diff --git a/include/gudhi/graph_simplicial_complex.h b/include/gudhi/graph_simplicial_complex.h index d84421b2..49fe56cc 100644 --- a/include/gudhi/graph_simplicial_complex.h +++ b/include/gudhi/graph_simplicial_complex.h @@ -4,7 +4,7 @@ * * Author(s): Clément Maria * - * Copyright (C) 2014 INRIA + * 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 @@ -42,6 +42,12 @@ struct vertex_filtration_t { typedef boost::vertex_property_tag kind; }; +/** \brief Proximity_graph contains the vertices and edges with their filtration values in order to store the result + * of `Gudhi::compute_proximity_graph` function. + * + * \tparam SimplicialComplexForProximityGraph furnishes `Filtration_value` type definition. + * + */ template <typename SimplicialComplexForProximityGraph> using Proximity_graph = typename boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS , boost::property < vertex_filtration_t, typename SimplicialComplexForProximityGraph::Filtration_value > diff --git a/include/gudhi/pick_n_random_points.h b/include/gudhi/pick_n_random_points.h index 8c90b6bf..64821e5d 100644 --- a/include/gudhi/pick_n_random_points.h +++ b/include/gudhi/pick_n_random_points.h @@ -4,7 +4,7 @@ * * Author(s): Siargey Kachanovich * - * Copyright (C) 2016 INRIA + * 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 diff --git a/include/gudhi/random_point_generators.h b/include/gudhi/random_point_generators.h index 9df77760..f8107c8b 100644 --- a/include/gudhi/random_point_generators.h +++ b/include/gudhi/random_point_generators.h @@ -4,7 +4,7 @@ * * Author(s): Clement Jamin * - * Copyright (C) 2016 INRIA + * 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 @@ -190,7 +190,8 @@ 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>()) { + 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) { diff --git a/include/gudhi/read_persistence_from_file.h b/include/gudhi/read_persistence_from_file.h index 83b89d0e..4a2b9d68 100644 --- a/include/gudhi/read_persistence_from_file.h +++ b/include/gudhi/read_persistence_from_file.h @@ -4,7 +4,7 @@ * * Author(s): Pawel Dlotko * - * Copyright (C) 2016 INRIA (France) + * 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 @@ -57,7 +57,7 @@ std::vector<std::pair<double, double> > read_persistence_intervals_in_one_dimens std::string line; std::vector<std::pair<double, double> > barcode_initial = - read_persistence_intervals_in_dimension(filename, (int)dimension); + read_persistence_intervals_in_dimension(filename, static_cast<int>(dimension)); std::vector<std::pair<double, double> > final_barcode; final_barcode.reserve(barcode_initial.size()); @@ -92,8 +92,8 @@ std::vector<std::pair<double, double> > read_persistence_intervals_in_one_dimens if ((barcode_initial[i].second == std::numeric_limits<double>::infinity()) && (what_to_substitute_for_infinite_bar != -1)) { - if (barcode_initial[i].first < what_to_substitute_for_infinite_bar) // if only birth < death. - { + if (barcode_initial[i].first < what_to_substitute_for_infinite_bar) { + // if only birth < death. final_barcode.push_back( std::pair<double, double>(barcode_initial[i].first, what_to_substitute_for_infinite_bar)); } diff --git a/include/gudhi/reader_utils.h b/include/gudhi/reader_utils.h index 90be4fc7..26eeb76d 100644 --- a/include/gudhi/reader_utils.h +++ b/include/gudhi/reader_utils.h @@ -4,7 +4,7 @@ * * Author(s): Clement Maria, Pawel Dlotko, Clement Jamin * - * Copyright (C) 2014 INRIA + * 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 diff --git a/include/gudhi/sparsify_point_set.h b/include/gudhi/sparsify_point_set.h index 7d3b97fb..db10e0b1 100644 --- a/include/gudhi/sparsify_point_set.h +++ b/include/gudhi/sparsify_point_set.h @@ -4,7 +4,7 @@ * * Author(s): Clement Jamin * - * Copyright (C) 2016 INRIA + * 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 diff --git a/include/gudhi/writing_persistence_to_file.h b/include/gudhi/writing_persistence_to_file.h new file mode 100644 index 00000000..34448576 --- /dev/null +++ b/include/gudhi/writing_persistence_to_file.h @@ -0,0 +1,117 @@ +/* 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): Pawel Dlotko + * + * Copyright (C) 2017 Swansea University, UK + * + * 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 WRITING_PERSISTENCE_TO_FILE_H_ +#define WRITING_PERSISTENCE_TO_FILE_H_ + +#include <iostream> +#include <string> +#include <limits> + +namespace Gudhi { + +/** +* This is a class to store persistence intervals. Its main purpose is to +* exchange data in between different packages and provide unified way +* of writing a collection of persistence intervals to file. +**/ +template <typename Filtration_type, typename Coefficient_field> +class Persistence_interval_common { + public: + /** + * Constructor taking as an input birth and death of the pair. + **/ + Persistence_interval_common(Filtration_type birth, Filtration_type death) + : birth_(birth), + death_(death), + dimension_(std::numeric_limits<unsigned>::max()), + arith_element_(std::numeric_limits<Coefficient_field>::max()) {} + + /** + * Constructor taking as an input birth, death and dimension of the pair. + **/ + Persistence_interval_common(Filtration_type birth, Filtration_type death, unsigned dim) + : birth_(birth), death_(death), dimension_(dim), arith_element_(std::numeric_limits<Coefficient_field>::max()) {} + + /** +* Constructor taking as an input birth, death, dimension of the pair as well +* as the number p such that this interval is present over Z_p field. +**/ + Persistence_interval_common(Filtration_type birth, Filtration_type death, unsigned dim, Coefficient_field field) + : birth_(birth), death_(death), dimension_(dim), arith_element_(field) {} + + /** + * Operator to compare two persistence pairs. During the comparision all the + * fields: birth, death, dimensiona and arith_element_ are taken into account + * and they all have to be equal for two pairs to be equal. + **/ + bool operator==(const Persistence_interval_common& i2) const { + return ((this->birth_ == i2.birth_) && (this->death_ == i2.death_) && (this->dimension_ == i2.dimension_) && + (this->arith_element_ == i2.arith_element_)); + } + + /** + * Check if two persistence paris are not equal. + **/ + bool operator!=(const Persistence_interval_common& i2) const { return (!((*this) == i2)); } + + /** + * Operator to compare objects of a type Persistence_interval_common. + * One intervals is smaller than the other if it has lower persistence. + * Note that this operator do not take Arith_element into account when doing comparisions. + **/ + bool operator<(const Persistence_interval_common& i2) const { + return fabs(this->death_ - this->birth_) < fabs(i2.death_ - i2.birth_); + } + + friend std::ostream& operator<<(std::ostream& out, const Persistence_interval_common& it) { + if (it.arith_element_ != std::numeric_limits<Coefficient_field>::max()) { + out << it.arith_element_ << " "; + } + if (it.dimension_ != std::numeric_limits<unsigned>::max()) { + out << it.dimension_ << " "; + } + out << it.birth_ << " " << it.death_ << " "; + return out; + } + + private: + Filtration_type birth_; + Filtration_type death_; + unsigned dimension_; + Coefficient_field arith_element_; +}; + +/** + * This function write a vector<Persistence_interval_common> to a stream +**/ +template <typename Persistence_interval_range> +void write_persistence_intervals_to_stream(const Persistence_interval_range& intervals, + std::ostream& out = std::cout) { + for (auto interval : intervals) { + out << interval << "\n"; + } +} + +} // namespace Gudhi + +#endif // WRITING_PERSISTENCE_TO_FILE_H_ |