diff options
| author | ROUVREAU Vincent <vincent.rouvreau@inria.fr> | 2019-06-12 07:46:27 +0200 |
|---|---|---|
| committer | ROUVREAU Vincent <vincent.rouvreau@inria.fr> | 2019-06-12 07:46:27 +0200 |
| commit | 6c57ce009e242db0554aaf9316750dd0bb693d6a (patch) | |
| tree | 94409a8ae609c48a6212ff61a90f7c083e0d8fcf /src | |
| parent | 4e33acbac51c8c2348dc88a16eb38e14c8ef724a (diff) | |
| parent | f58f0bb2cb99076d0cd3a11ad39f3277213e3f5e (diff) | |
Merge branch 'master' into main_documentation_pages
Diffstat (limited to 'src')
31 files changed, 964 insertions, 278 deletions
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h index 32dfcc16..0bf12b1a 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h @@ -58,7 +58,7 @@ #include <limits> // for numeric_limits<> #if CGAL_VERSION_NR < 1041101000 -// Make compilation fail - required for external projects - https://gitlab.inria.fr/GUDHI/gudhi-devel/issues/10 +// Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10 # error Alpha_complex_3d is only available for CGAL >= 4.11 #endif diff --git a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp index 622fcae8..b46b6da5 100644 --- a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp +++ b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp @@ -42,7 +42,7 @@ // Use dynamic_dimension_tag for the user to be able to set dimension typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kernel_d; // Use static dimension_tag for the user not to be able to set dimension -typedef CGAL::Epick_d< CGAL::Dimension_tag<2> > Kernel_s; +typedef CGAL::Epick_d< CGAL::Dimension_tag<3> > Kernel_s; // The triangulation uses the default instantiation of the TriangulationDataStructure template parameter typedef boost::mpl::list<Kernel_d, Kernel_s> list_of_kernel_variants; diff --git a/src/Alpha_complex/utilities/alphacomplex.md b/src/Alpha_complex/utilities/alphacomplex.md index 50a39d32..fcd16a3b 100644 --- a/src/Alpha_complex/utilities/alphacomplex.md +++ b/src/Alpha_complex/utilities/alphacomplex.md @@ -33,7 +33,7 @@ a prime number). where
`<input OFF file>` is the path to the input point cloud in
-[nOFF ASCII format](http://www.geomview.org/docs/html/OFF.html).
+[nOFF ASCII format]({{ site.officialurl }}/doc/latest/fileformats.html#FileFormatsOFF).
**Allowed options**
@@ -87,7 +87,7 @@ a prime number). ```
where `<input OFF file>` is the path to the input point cloud in
-[nOFF ASCII format](http://www.geomview.org/docs/html/OFF.html).
+[nOFF ASCII format]({{ site.officialurl }}/doc/latest/fileformats.html#FileFormatsOFF).
**Allowed options**
@@ -103,8 +103,10 @@ to be recorded. Enter a negative value to see zero length intervals. * `-c [ --cuboid-file ]` is the path to the file describing the periodic domain.
It must be in the format described
[here]({{ site.officialurl }}/doc/latest/fileformats.html#FileFormatsIsoCuboid).
+Default version is not periodic.
* `-w [ --weight-file ]` is the path to the file containing the weights of the
points (one value per line).
+Default version is not weighted.
* `-e [ --exact ]` for the exact computation version (not compatible with
weight and periodic version).
* `-f [ --fast ]` for the fast computation version.
diff --git a/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h b/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h index f8fce96c..6fd058a8 100644 --- a/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h +++ b/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h @@ -41,6 +41,9 @@ namespace persistence_diagram { * * \image html perturb_pd.png On this picture, the red edges represent the matching. The bottleneck distance is the length of the longest edge. * + * This implementation is based on ideas from "Geometry Helps in Bottleneck Matching and Related Problems" + * \cite DBLP:journals/algorithmica/EfratIK01. Another relevant publication, although it was not used is + * "Geometry Helps to Compare Persistence Diagrams" \cite Kerber:2017:GHC:3047249.3064175. */ /** @} */ // end defgroup bottleneck_distance diff --git a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h index 36a63ea0..8c12d353 100644 --- a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h +++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h @@ -33,6 +33,7 @@ #include <unordered_set> #include <vector> #include <algorithm> // for std::max +#include <cmath> // for std::abs namespace Gudhi { @@ -45,7 +46,7 @@ struct Square_query { typedef Internal_point Point_d; typedef double FT; bool contains(Point_d p) const { - return std::max(std::abs(p.x()-c.x()), std::abs(p.y()-c.y())) <= size; + return (std::max)(std::abs(p.x()-c.x()), std::abs(p.y()-c.y())) <= size; } bool inner_range_intersects(CGAL::Kd_tree_rectangle<FT, D> const&r) const { return diff --git a/src/Doxyfile.in b/src/Doxyfile.in index 54a438d4..f80d4505 100644 --- a/src/Doxyfile.in +++ b/src/Doxyfile.in @@ -821,7 +821,8 @@ EXCLUDE_SYMBOLS = EXAMPLE_PATH = biblio/ \ example/ \ - utilities/ + utilities/ \ + data/ # If the value of the EXAMPLE_PATH tag contains directories, you can use the # EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp and @@ -1440,7 +1441,7 @@ MATHJAX_FORMAT = HTML-CSS # The default value is: http://cdn.mathjax.org/mathjax/latest. # This tag requires that the tag USE_MATHJAX is set to YES. -MATHJAX_RELPATH = ../common +MATHJAX_RELPATH = https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.2 # The MATHJAX_EXTENSIONS tag can be used to specify one or more MathJax # extension names that should be enabled during MathJax rendering. For example diff --git a/src/Nerve_GIC/doc/Intro_graph_induced_complex.h b/src/Nerve_GIC/doc/Intro_graph_induced_complex.h index bc8aecc3..e72d63dd 100644 --- a/src/Nerve_GIC/doc/Intro_graph_induced_complex.h +++ b/src/Nerve_GIC/doc/Intro_graph_induced_complex.h @@ -37,8 +37,7 @@ namespace cover_complex { * neato (from <a target="_blank" href="http://www.graphviz.org/">graphviz</a>), * <a target="_blank" href="http://www.geomview.org/">geomview</a>, * <a target="_blank" href="https://github.com/MLWave/kepler-mapper">KeplerMapper</a>. - * Input point clouds are assumed to be - * <a target="_blank" href="http://www.geomview.org/docs/html/OFF.html">OFF files</a>. + * Input point clouds are assumed to be \ref FileFormatsOFF "OFF files" * * \section covers Covers * diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h index c3085dff..11bb4e85 100644 --- a/src/Nerve_GIC/include/gudhi/GIC.h +++ b/src/Nerve_GIC/include/gudhi/GIC.h @@ -52,7 +52,7 @@ #include <string> #include <limits> // for numeric_limits #include <utility> // for std::pair<> -#include <algorithm> // for std::max +#include <algorithm> // for (std::max) #include <random> #include <cassert> #include <cmath> @@ -457,7 +457,7 @@ class Cover_complex { template <typename Distance> 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); + m = (std::min)(m, n - 1); double delta = 0; if (verbose) std::cout << n << " points in R^" << data_dimension << std::endl; @@ -475,8 +475,8 @@ class Cover_complex { 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); + 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; @@ -489,8 +489,8 @@ class Cover_complex { 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); + 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; } @@ -586,7 +586,7 @@ class Cover_complex { if (type == "GIC") { boost::graph_traits<Graph>::edge_iterator ei, ei_end; for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei) - reso = std::max(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] - + reso = (std::max)(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] - func[index[boost::target(*ei, one_skeleton)]])); if (verbose) std::cout << "resolution = " << reso << std::endl; resolution_double = reso; @@ -595,7 +595,7 @@ class Cover_complex { if (type == "Nerve") { boost::graph_traits<Graph>::edge_iterator ei, ei_end; for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei) - reso = std::max(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] - + reso = (std::max)(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] - func[index[boost::target(*ei, one_skeleton)]]) / gain); if (verbose) std::cout << "resolution = " << reso << std::endl; @@ -640,11 +640,11 @@ class Cover_complex { } // Read function values and compute min and max - double minf = std::numeric_limits<float>::max(); + double minf = (std::numeric_limits<float>::max)(); double maxf = std::numeric_limits<float>::lowest(); for (int i = 0; i < n; i++) { - minf = std::min(minf, func[i]); - maxf = std::max(maxf, func[i]); + minf = (std::min)(minf, func[i]); + maxf = (std::max)(maxf, func[i]); } if (verbose) std::cout << "Min function value = " << minf << " and Max function value = " << maxf << std::endl; @@ -899,7 +899,7 @@ class Cover_complex { Weight_map weight = boost::get(boost::edge_weight, one_skeleton); Index_map index = boost::get(boost::vertex_index, one_skeleton); std::vector<double> mindist(n); - for (int j = 0; j < n; j++) mindist[j] = std::numeric_limits<double>::max(); + for (int j = 0; j < n; j++) mindist[j] = (std::numeric_limits<double>::max)(); // Compute the geodesic distances to subsamples with Dijkstra #ifdef GUDHI_USE_TBB @@ -1027,10 +1027,10 @@ class Cover_complex { std::ofstream graphic(mapp); double maxv = std::numeric_limits<double>::lowest(); - double minv = std::numeric_limits<double>::max(); + double minv = (std::numeric_limits<double>::max)(); for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) { - maxv = std::max(maxv, iit->second.second); - minv = std::min(minv, iit->second.second); + maxv = (std::max)(maxv, iit->second.second); + minv = (std::min)(minv, iit->second.second); } int k = 0; @@ -1162,10 +1162,10 @@ class Cover_complex { // Compute max and min double maxf = std::numeric_limits<double>::lowest(); - double minf = std::numeric_limits<double>::max(); + double minf = (std::numeric_limits<double>::max)(); for (std::map<int, double>::iterator it = cover_std.begin(); it != cover_std.end(); it++) { - maxf = std::max(maxf, it->second); - minf = std::min(minf, it->second); + maxf = (std::max)(maxf, it->second); + minf = (std::min)(minf, it->second); } // Build filtration @@ -1294,8 +1294,8 @@ class Cover_complex { * */ double compute_p_value() { - double distancemin = std::numeric_limits<double>::max(); int N = PD.size(); - for (int i = 0; i < N; i++) distancemin = std::min(distancemin, 0.5 * std::abs(PD[i].second - PD[i].first)); + double distancemin = (std::numeric_limits<double>::max)(); int N = PD.size(); + for (int i = 0; i < N; i++) distancemin = (std::min)(distancemin, 0.5 * std::abs(PD[i].second - PD[i].first)); double p_value = 1 - compute_confidence_level_from_distance(distancemin); if (verbose) std::cout << "p value = " << p_value << std::endl; return p_value; @@ -1317,7 +1317,7 @@ class Cover_complex { for (auto const& simplex : simplices) { int numvert = simplex.size(); double filt = std::numeric_limits<double>::lowest(); - for (int i = 0; i < numvert; i++) filt = std::max(cover_color[simplex[i]].second, filt); + for (int i = 0; i < numvert; i++) filt = (std::max)(cover_color[simplex[i]].second, filt); complex.insert_simplex_and_subfaces(simplex, filt); if (dimension < simplex.size() - 1) dimension = simplex.size() - 1; } @@ -1361,8 +1361,8 @@ class Cover_complex { int vt = cover[index[boost::target(*ei, one_skeleton)]][j]; if (cover_fct[vs] == cover_fct[vt] + 1 || cover_fct[vt] == cover_fct[vs] + 1) { std::vector<int> edge(2); - edge[0] = std::min(vs, vt); - edge[1] = std::max(vs, vt); + edge[0] = (std::min)(vs, vt); + edge[1] = (std::max)(vs, vt); simplices.push_back(edge); goto afterLoop; } diff --git a/src/Persistence_representations/doc/Persistence_representations_doc.h b/src/Persistence_representations/doc/Persistence_representations_doc.h index 4d850a02..668904c9 100644 --- a/src/Persistence_representations/doc/Persistence_representations_doc.h +++ b/src/Persistence_representations/doc/Persistence_representations_doc.h @@ -227,6 +227,19 @@ namespace Persistence_representations { to diagonal are given then sometimes the kernel have support that reaches the region below the diagonal. If the value of this parameter is true, then the values below diagonal can be erased. + In addition to the previous method, we also provide two more methods to perform exact calculations, in the sense that we use functions + instead of matrices to define the kernel between the points of the diagrams. + Indeed, in both of these exact methods, the kernel is no longer provided as a square matrix, or a filter (see parameters above), but rather as + a function assigning a real value to a pair of points in the plane. + + In the first of these exact methods, we aim at obtaining a finite-dimensional representation of the diagram, so we still use a grid of pixels. + On the other hand, in the second exact method, we represent diagrams implicitly as functions (i.e. infinite-dimensional representations). This way, we no longer require grids, + but only scalar products and distances are available with these implicit representations. This type of representations is known as + kernel methods (see \ref sec_persistence_kernels below for more details on kernels). + + Names can be a bit confusing so we recall that, with this second exact method, we implicitly define a kernel representation of diagrams that is built from a kernel between points + in the plane. Hence, we have two kernels here, which are independent. One is defined between points in the plane (its type in the code is Kernel2D), and is a template parameter, + whereas the other is defined between persistence diagrams (it is the scalar product of the infinite-dimensional representations of the diagrams). \section sec_persistence_vectors Persistence vectors <b>Reference manual:</b> \ref Gudhi::Persistence_representations::Vector_distances_in_diagram <br> @@ -250,6 +263,35 @@ namespace Persistence_representations { absolute value of differences between coordinates. A scalar product is a sum of products of values at the corresponding positions of two vectors. + \section sec_persistence_kernels Kernels on persistence diagrams + <b>Reference manual:</b> \ref Gudhi::Persistence_representations::Sliced_Wasserstein <br> + <b>Reference manual:</b> \ref Gudhi::Persistence_representations::Persistence_heat_maps <br> + + Kernels for persistence diagrams can be regarded as infinite-dimensional vectorizations. More specifically, + they are similarity functions whose evaluations on pairs of persistence diagrams equals the scalar products + between images of these pairs under a map \f$\Phi\f$ taking values in a specific (possibly non Euclidean) Hilbert space \f$k(D_i, D_j) = \langle \Phi(D_i),\Phi(D_j)\rangle\f$. + Reciprocally, classical results of learning theory ensure that such a \f$\Phi\f$ exists for a given similarity function \f$k\f$ if and only if \f$k\f$ is <i>positive semi-definite</i>. + Kernels are designed for algorithms that can be <i>kernelized</i>, i.e., algorithms that only require to know scalar products between instances in order to run. + Examples of such algorithms include Support Vector Machines, Principal Component Analysis and Ridge Regression. + + There have been several attempts at defining kernels, i.e., positive semi-definite functions, between persistence diagrams within the last few years. We provide implementation + for the <i>Sliced Wasserstein kernel</i>---see \cite pmlr-v70-carriere17a, which takes the form of a Gaussian kernel with a specific distance between persistence diagrams + called the <i>Sliced Wasserstein distance</i>: \f$k(D_1,D_2)={\rm exp}\left(-\frac{SW(D_1,D_2)}{2\sigma^2}\right)\f$. Other kernels such as the Persistence Weighted Gaussian kernel or + the Persistence Scale Space kernel are implemented in Persistence_heat_maps. + + When launching: + + \code $> ./Sliced_Wasserstein + \endcode + + the program output is: + + \code $> Approx SW kernel: 0.0693743 + $> Exact SW kernel: 0.0693218 + $> Distance induced by approx SW kernel: 1.36428 + $> Distance induced by exact SW kernel: 1.3643 + \endcode + */ /** @} */ // end defgroup Persistence_representations diff --git a/src/Persistence_representations/example/CMakeLists.txt b/src/Persistence_representations/example/CMakeLists.txt index 33558df3..a7c6ef39 100644 --- a/src/Persistence_representations/example/CMakeLists.txt +++ b/src/Persistence_representations/example/CMakeLists.txt @@ -26,3 +26,7 @@ add_test(NAME Persistence_representations_example_heat_maps COMMAND $<TARGET_FILE:Persistence_representations_example_heat_maps>) install(TARGETS Persistence_representations_example_heat_maps DESTINATION bin) +add_executable ( Sliced_Wasserstein sliced_wasserstein.cpp ) +add_test(NAME Sliced_Wasserstein + COMMAND $<TARGET_FILE:Sliced_Wasserstein>) +install(TARGETS Sliced_Wasserstein DESTINATION bin) diff --git a/src/Persistence_representations/example/persistence_heat_maps.cpp b/src/Persistence_representations/example/persistence_heat_maps.cpp index 323b57e9..45208b68 100644 --- a/src/Persistence_representations/example/persistence_heat_maps.cpp +++ b/src/Persistence_representations/example/persistence_heat_maps.cpp @@ -2,9 +2,12 @@ * (Geometric Understanding in Higher Dimensions) is a generic C++ * library for computational topology. * - * Author(s): Pawel Dlotko + * Author(s): Pawel Dlotko and Mathieu Carriere * - * Copyright (C) 2016 Inria + * Copyright (C) 2019 Inria + * + * Modifications: + * - 2018/04 MC: Add persistence heat maps computation * * 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 @@ -21,10 +24,20 @@ */ #include <gudhi/Persistence_heat_maps.h> +#include <gudhi/common_persistence_representations.h> #include <iostream> #include <vector> #include <utility> +#include <functional> +#include <cmath> + +std::function<double(std::pair<double, double>, std::pair<double, double>)> Gaussian_function(double sigma) { + return [=](std::pair<double, double> p, std::pair<double, double> q) { + return std::exp(-((p.first - q.first) * (p.first - q.first) + (p.second - q.second) * (p.second - q.second)) / + (sigma)); + }; +} using constant_scaling_function = Gudhi::Persistence_representations::constant_scaling_function; using Persistence_heat_maps = Gudhi::Persistence_representations::Persistence_heat_maps<constant_scaling_function>; @@ -76,5 +89,13 @@ int main(int argc, char** argv) { // to compute scalar product of hm1 and hm2: std::cout << "Scalar product is : " << hm1.compute_scalar_product(hm2) << std::endl; + Persistence_heat_maps hm1k(persistence1, Gaussian_function(1.0)); + Persistence_heat_maps hm2k(persistence2, Gaussian_function(1.0)); + Persistence_heat_maps hm1i(persistence1, Gaussian_function(1.0), 20, 20, 0, 11, 0, 11); + Persistence_heat_maps hm2i(persistence2, Gaussian_function(1.0), 20, 20, 0, 11, 0, 11); + std::cout << "Scalar product computed with exact 2D kernel on grid is : " << hm1i.compute_scalar_product(hm2i) + << std::endl; + std::cout << "Scalar product computed with exact 2D kernel is : " << hm1k.compute_scalar_product(hm2k) << std::endl; + return 0; } diff --git a/src/Persistence_representations/example/sliced_wasserstein.cpp b/src/Persistence_representations/example/sliced_wasserstein.cpp new file mode 100644 index 00000000..089172a0 --- /dev/null +++ b/src/Persistence_representations/example/sliced_wasserstein.cpp @@ -0,0 +1,59 @@ +/* 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): Mathieu Carriere + * + * Copyright (C) 2018 INRIA (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <gudhi/Sliced_Wasserstein.h> + +#include <iostream> +#include <vector> +#include <utility> + +using Persistence_diagram = Gudhi::Persistence_representations::Persistence_diagram; +using SW = Gudhi::Persistence_representations::Sliced_Wasserstein; + +int main(int argc, char** argv) { + + Persistence_diagram persistence1, persistence2; + + persistence1.push_back(std::make_pair(1, 2)); + persistence1.push_back(std::make_pair(6, 8)); + persistence1.push_back(std::make_pair(0, 4)); + persistence1.push_back(std::make_pair(3, 8)); + + persistence2.push_back(std::make_pair(2, 9)); + persistence2.push_back(std::make_pair(1, 6)); + persistence2.push_back(std::make_pair(3, 5)); + persistence2.push_back(std::make_pair(6, 10)); + + + SW sw1(persistence1, 1, 100); + SW sw2(persistence2, 1, 100); + + SW swex1(persistence1, 1, -1); + SW swex2(persistence2, 1, -1); + + std::cout << "Approx SW kernel: " << sw1.compute_scalar_product(sw2) << std::endl; + std::cout << "Exact SW kernel: " << swex1.compute_scalar_product(swex2) << std::endl; + std::cout << "Distance induced by approx SW kernel: " << sw1.distance(sw2) << std::endl; + std::cout << "Distance induced by exact SW kernel: " << swex1.distance(swex2) << std::endl; + + return 0; +} diff --git a/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h b/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h index 35e51e63..a8458bda 100644 --- a/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h +++ b/src/Persistence_representations/include/gudhi/Persistence_heat_maps.h @@ -2,9 +2,12 @@ * (Geometric Understanding in Higher Dimensions) is a generic C++ * library for computational topology. * - * Author(s): Pawel Dlotko + * Author(s): Pawel Dlotko and Mathieu Carriere * - * Copyright (C) 2016 Inria + * Modifications: + * - 2018/04 MC: Add discrete/non-discrete mechanism and non-discrete version + * + * Copyright (C) 2019 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 @@ -44,7 +47,7 @@ namespace Persistence_representations { /** * This is a simple procedure to create n by n (or 2*pixel_radius times 2*pixel_radius cubical approximation of a *Gaussian kernel. -**/ + **/ std::vector<std::vector<double> > create_Gaussian_filter(size_t pixel_radius, double sigma) { bool dbg = false; // we are computing the kernel mask to 2 standard deviations away from the center. We discretize it in a grid of a @@ -74,7 +77,7 @@ std::vector<std::vector<double> > create_Gaussian_filter(size_t pixel_radius, do for (int y = -pixel_radius; y <= static_cast<int>(pixel_radius); y++) { double real_x = 2 * sigma * x / pixel_radius; double real_y = 2 * sigma * y / pixel_radius; - r = sqrt(real_x * real_x + real_y * real_y); + r = std::sqrt(real_x * real_x + real_y * real_y); kernel[x + pixel_radius][y + pixel_radius] = (exp(-(r * r) / sigma_sqr)) / (3.141592 * sigma_sqr); sum += kernel[x + pixel_radius][y + pixel_radius]; } @@ -100,18 +103,18 @@ std::vector<std::vector<double> > create_Gaussian_filter(size_t pixel_radius, do } /* -* There are various options to scale the points depending on their location. One can for instance: -* (1) do nothing (scale all of them with the weight 1), as in the function constant_function -* (2) Scale them by the distance to the diagonal. This is implemented in function -* (3) Scale them with the square of their distance to diagonal. This is implemented in function -* (4) Scale them with -*/ + * There are various options to scale the points depending on their location. One can for instance: + * (1) do nothing (scale all of them with the weight 1), as in the function constant_function + * (2) Scale them by the distance to the diagonal. This is implemented in function + * (3) Scale them with the square of their distance to diagonal. This is implemented in function + * (4) Scale them with + */ /** * This is one of a scaling functions used to weight points depending on their persistence and/or location in the *diagram. * This particular functionality is a function which always assign value 1 to a point in the diagram. -**/ + **/ class constant_scaling_function { public: double operator()(const std::pair<double, double>& point_in_diagram) { return 1; } @@ -121,13 +124,13 @@ class constant_scaling_function { * This is one of a scaling functions used to weight points depending on their persistence and/or location in the *diagram. * The scaling given by this function to a point (b,d) is Euclidean distance of (b,d) from diagonal. -**/ + **/ class distance_from_diagonal_scaling { public: double operator()(const std::pair<double, double>& point_in_diagram) { // (point_in_diagram.first+point_in_diagram.second)/2.0 - return sqrt(pow((point_in_diagram.first - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2) + - pow((point_in_diagram.second - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2)); + return std::sqrt(std::pow((point_in_diagram.first - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2) + + std::pow((point_in_diagram.second - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2)); } }; @@ -135,12 +138,12 @@ class distance_from_diagonal_scaling { * This is one of a scaling functions used to weight points depending on their persistence and/or location in the *diagram. * The scaling given by this function to a point (b,d) is a square of Euclidean distance of (b,d) from diagonal. -**/ + **/ class squared_distance_from_diagonal_scaling { public: double operator()(const std::pair<double, double>& point_in_diagram) { - return pow((point_in_diagram.first - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2) + - pow((point_in_diagram.second - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2); + return std::pow((point_in_diagram.first - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2) + + std::pow((point_in_diagram.second - (point_in_diagram.first + point_in_diagram.second) / 2.0), 2); } }; @@ -148,7 +151,7 @@ class squared_distance_from_diagonal_scaling { * This is one of a scaling functions used to weight points depending on their persistence and/or location in the *diagram. * The scaling given by this function to a point (b,d) is an arctan of a persistence of a point (i.e. arctan( b-d ). -**/ + **/ class arc_tan_of_persistence_of_point { public: double operator()(const std::pair<double, double>& point_in_diagram) { @@ -162,16 +165,16 @@ class arc_tan_of_persistence_of_point { * This scaling function do not only depend on a point (p,d) in the diagram, but it depends on the whole diagram. * The longest persistence pair get a scaling 1. Any other pair get a scaling belong to [0,1], which is proportional * to the persistence of that pair. -**/ + **/ class weight_by_setting_maximal_interval_to_have_length_one { public: - weight_by_setting_maximal_interval_to_have_length_one(double len) : letngth_of_maximal_interval(len) {} + weight_by_setting_maximal_interval_to_have_length_one(double len) : length_of_maximal_interval(len) {} double operator()(const std::pair<double, double>& point_in_diagram) { - return (point_in_diagram.second - point_in_diagram.first) / this->letngth_of_maximal_interval; + return (point_in_diagram.second - point_in_diagram.first) / this->length_of_maximal_interval; } private: - double letngth_of_maximal_interval; + double length_of_maximal_interval; }; /** @@ -179,7 +182,7 @@ class weight_by_setting_maximal_interval_to_have_length_one { * \brief A class implementing persistence heat maps. * * \ingroup Persistence_representations -**/ + **/ // This class implements the following concepts: Vectorized_topological_data, Topological_data_with_distances, // Real_valued_topological_data, Topological_data_with_averages, Topological_data_with_scalar_product @@ -189,7 +192,7 @@ class Persistence_heat_maps { /** * The default constructor. A scaling function from the diagonal is set up to a constant function. The image is not *erased below the diagonal. The Gaussian have diameter 5. - **/ + **/ Persistence_heat_maps() { Scalling_of_kernels f; this->f = f; @@ -210,7 +213,7 @@ class Persistence_heat_maps { *std::numeric_limits<double>::max(), in which case the program compute the values based on the data, * (6) a max x and y value of points that are to be taken into account. By default it is set to *std::numeric_limits<double>::max(), in which case the program compute the values based on the data. - **/ + **/ Persistence_heat_maps(const std::vector<std::pair<double, double> >& interval, std::vector<std::vector<double> > filter = create_Gaussian_filter(5, 1), bool erase_below_diagonal = false, size_t number_of_pixels = 1000, @@ -218,12 +221,12 @@ class Persistence_heat_maps { double max_ = std::numeric_limits<double>::max()); /** - * Construction that takes at the input a name of a file with persistence intervals, a filter (radius 5 by - *default), a scaling function (constant by default), a boolean value which determines if the area of image below - *diagonal should, or should not be erased (should by default). The next parameter is the number of pixels in each - *direction (set to 1000 by default) and min and max values of images (both set to std::numeric_limits<double>::max() - *by default. If this is the case, the program will pick the right values based on the data). - **/ + * Construction that takes at the input a name of a file with persistence intervals, a filter (radius 5 by + *default), a scaling function (constant by default), a boolean value which determines if the area of image below + *diagonal should, or should not be erased (should by default). The next parameter is the number of pixels in each + *direction (set to 1000 by default) and min and max values of images (both set to std::numeric_limits<double>::max() + *by default. If this is the case, the program will pick the right values based on the data). + **/ /** * Construction that takes at the input the following parameters: * (1) A name of a file with persistence intervals. The file should be readable by the function @@ -237,7 +240,7 @@ class Persistence_heat_maps { *std::numeric_limits<double>::max(), in which case the program compute the values based on the data, * (6) a max x and y value of points that are to be taken into account. By default it is set to *std::numeric_limits<double>::max(), in which case the program compute the values based on the data. - **/ + **/ Persistence_heat_maps(const char* filename, std::vector<std::vector<double> > filter = create_Gaussian_filter(5, 1), bool erase_below_diagonal = false, size_t number_of_pixels = 1000, double min_ = std::numeric_limits<double>::max(), @@ -245,22 +248,37 @@ class Persistence_heat_maps { unsigned dimension = std::numeric_limits<unsigned>::max()); /** + * Construction that takes as inputs (1) the diagram, and (2) the grid parameters (min, max and number of samples for + * x and y axes) + **/ + Persistence_heat_maps(const std::vector<std::pair<double, double> >& interval, + const std::function<double(std::pair<double, double>, std::pair<double, double>)>& kernel, + size_t number_of_x_pixels, size_t number_of_y_pixels, double min_x = 0, double max_x = 1, + double min_y = 0, double max_y = 1); + + /** + * Construction that takes only the diagram as input (weight and 2D kernel are template parameters) + **/ + Persistence_heat_maps(const std::vector<std::pair<double, double> >& interval, + const std::function<double(std::pair<double, double>, std::pair<double, double>)>& kernel); + + /** * Compute a mean value of a collection of heat maps and store it in the current object. Note that all the persistence *maps send in a vector to this procedure need to have the same parameters. * If this is not the case, the program will throw an exception. - **/ + **/ void compute_mean(const std::vector<Persistence_heat_maps*>& maps); /** * Compute a median value of a collection of heat maps and store it in the current object. Note that all the *persistence maps send in a vector to this procedure need to have the same parameters. * If this is not the case, the program will throw an exception. - **/ + **/ void compute_median(const std::vector<Persistence_heat_maps*>& maps); /** * Compute a percentage of active (i.e) values above the cutoff of a collection of heat maps. - **/ + **/ void compute_percentage_of_active(const std::vector<Persistence_heat_maps*>& maps, size_t cutoff = 1); // put to file subroutine @@ -268,18 +286,18 @@ class Persistence_heat_maps { * The function outputs the persistence image to a text file. The format as follow: * In the first line, the values min and max of the image are stored * In the next lines, we have the persistence images in a form of a bitmap image. - **/ + **/ void print_to_file(const char* filename) const; /** * A function that load a heat map from file to the current object (and erase whatever was stored in the current *object before). - **/ + **/ void load_from_file(const char* filename); /** * The procedure checks if min_, max_ and this->heat_maps sizes are the same. - **/ + **/ inline bool check_if_the_same(const Persistence_heat_maps& second) const { bool dbg = false; if (this->heat_map.size() != second.heat_map.size()) { @@ -302,17 +320,17 @@ class Persistence_heat_maps { /** * Return minimal range value of persistent image. - **/ + **/ inline double get_min() const { return this->min_; } /** * Return maximal range value of persistent image. - **/ + **/ inline double get_max() const { return this->max_; } /** * Operator == to check if to persistence heat maps are the same. - **/ + **/ bool operator==(const Persistence_heat_maps& rhs) const { bool dbg = false; if (!this->check_if_the_same(rhs)) { @@ -335,12 +353,12 @@ class Persistence_heat_maps { /** * Operator != to check if to persistence heat maps are different. - **/ + **/ bool operator!=(const Persistence_heat_maps& rhs) const { return !((*this) == rhs); } /** * A function to generate a gnuplot script to visualize the persistent image. - **/ + **/ void plot(const char* filename) const; template <typename Operation_type> @@ -370,7 +388,7 @@ class Persistence_heat_maps { /** * Multiplication of Persistence_heat_maps by scalar (so that all values of the heat map gets multiplied by that *scalar). - **/ + **/ Persistence_heat_maps multiply_by_scalar(double scalar) const { Persistence_heat_maps result; result.min_ = this->min_; @@ -389,56 +407,56 @@ class Persistence_heat_maps { /** * This function computes a sum of two objects of a type Persistence_heat_maps. - **/ + **/ friend Persistence_heat_maps operator+(const Persistence_heat_maps& first, const Persistence_heat_maps& second) { return operation_on_pair_of_heat_maps(first, second, std::plus<double>()); } /** -* This function computes a difference of two objects of a type Persistence_heat_maps. -**/ + * This function computes a difference of two objects of a type Persistence_heat_maps. + **/ friend Persistence_heat_maps operator-(const Persistence_heat_maps& first, const Persistence_heat_maps& second) { return operation_on_pair_of_heat_maps(first, second, std::minus<double>()); } /** -* This function computes a product of an object of a type Persistence_heat_maps with real number. -**/ + * This function computes a product of an object of a type Persistence_heat_maps with real number. + **/ friend Persistence_heat_maps operator*(double scalar, const Persistence_heat_maps& A) { return A.multiply_by_scalar(scalar); } /** -* This function computes a product of an object of a type Persistence_heat_maps with real number. -**/ + * This function computes a product of an object of a type Persistence_heat_maps with real number. + **/ friend Persistence_heat_maps operator*(const Persistence_heat_maps& A, double scalar) { return A.multiply_by_scalar(scalar); } /** -* This function computes a product of an object of a type Persistence_heat_maps with real number. -**/ + * This function computes a product of an object of a type Persistence_heat_maps with real number. + **/ Persistence_heat_maps operator*(double scalar) { return this->multiply_by_scalar(scalar); } /** * += operator for Persistence_heat_maps. - **/ + **/ Persistence_heat_maps operator+=(const Persistence_heat_maps& rhs) { *this = *this + rhs; return *this; } /** - * -= operator for Persistence_heat_maps. - **/ + * -= operator for Persistence_heat_maps. + **/ Persistence_heat_maps operator-=(const Persistence_heat_maps& rhs) { *this = *this - rhs; return *this; } /** - * *= operator for Persistence_heat_maps. - **/ + * *= operator for Persistence_heat_maps. + **/ Persistence_heat_maps operator*=(double x) { *this = *this * x; return *this; } /** - * /= operator for Persistence_heat_maps. - **/ + * /= operator for Persistence_heat_maps. + **/ Persistence_heat_maps operator/=(double x) { if (x == 0) throw("In operator /=, division by 0. Program terminated."); *this = *this * (1 / x); @@ -448,14 +466,14 @@ class Persistence_heat_maps { // Implementations of functions for various concepts. /** - * This function produce a vector of doubles based on a persistence heat map. It is required in a concept + * This function produce a vector of doubles based on a persistence heat map. It is required in a concept * Vectorized_topological_data - */ + */ std::vector<double> vectorize(int number_of_function) const; /** - * This function return the number of functions that allows vectorization of persistence heat map. It is required - *in a concept Vectorized_topological_data. - **/ + * This function return the number of functions that allows vectorization of persistence heat map. It is required + *in a concept Vectorized_topological_data. + **/ size_t number_of_vectorize_functions() const { return this->number_of_functions_for_vectorization; } /** @@ -464,45 +482,45 @@ class Persistence_heat_maps { * At the moment this function is not tested, since it is quite likely to be changed in the future. Given this, when *using it, keep in mind that it * will be most likely changed in the next versions. - **/ + **/ double project_to_R(int number_of_function) const; /** * The function gives the number of possible projections to R. This function is required by the *Real_valued_topological_data concept. - **/ + **/ size_t number_of_projections_to_R() const { return this->number_of_functions_for_projections_to_reals; } /** - * A function to compute distance between persistence heat maps. - * The parameter of this function is a const reference to an object of a class Persistence_heat_maps. - * This function is required in Topological_data_with_distances concept. -* For max norm distance, set power to std::numeric_limits<double>::max() - **/ + * A function to compute distance between persistence heat maps. + * The parameter of this function is a const reference to an object of a class Persistence_heat_maps. + * This function is required in Topological_data_with_distances concept. + * For max norm distance, set power to std::numeric_limits<double>::max() + **/ double distance(const Persistence_heat_maps& second_, double power = 1) const; /** * A function to compute averaged persistence heat map, based on vector of persistence heat maps. * This function is required by Topological_data_with_averages concept. - **/ + **/ void compute_average(const std::vector<Persistence_heat_maps*>& to_average); /** - * A function to compute scalar product of persistence heat maps. - * The parameter of this function is a const reference to an object of a class Persistence_heat_maps. - * This function is required in Topological_data_with_scalar_product concept. - **/ + * A function to compute scalar product of persistence heat maps. + * The parameter of this function is a const reference to an object of a class Persistence_heat_maps. + * This function is required in Topological_data_with_scalar_product concept. + **/ double compute_scalar_product(const Persistence_heat_maps& second_) const; // end of implementation of functions needed for concepts. /** * The x-range of the persistence heat map. - **/ + **/ std::pair<double, double> get_x_range() const { return std::make_pair(this->min_, this->max_); } /** * The y-range of the persistence heat map. - **/ + **/ std::pair<double, double> get_y_range() const { return this->get_x_range(); } protected: @@ -512,7 +530,6 @@ class Persistence_heat_maps { size_t number_of_functions_for_projections_to_reals; void construct(const std::vector<std::pair<double, double> >& intervals_, std::vector<std::vector<double> > filter = create_Gaussian_filter(5, 1), - bool erase_below_diagonal = false, size_t number_of_pixels = 1000, double min_ = std::numeric_limits<double>::max(), double max_ = std::numeric_limits<double>::max()); @@ -521,6 +538,12 @@ class Persistence_heat_maps { this->number_of_functions_for_projections_to_reals = 1; } + // Boolean indicating if we are computing persistence image (true) or persistence weighted gaussian kernel (false) + bool discrete = true; + std::function<double(std::pair<double, double>, std::pair<double, double>)> kernel; + std::vector<std::pair<double, double> > interval; + std::vector<double> weights; + // data Scalling_of_kernels f; bool erase_below_diagonal; @@ -529,6 +552,45 @@ class Persistence_heat_maps { std::vector<std::vector<double> > heat_map; }; +template <typename Scalling_of_kernels> +Persistence_heat_maps<Scalling_of_kernels>::Persistence_heat_maps( + const std::vector<std::pair<double, double> >& interval, + const std::function<double(std::pair<double, double>, std::pair<double, double>)>& kernel, + size_t number_of_x_pixels, size_t number_of_y_pixels, double min_x, double max_x, double min_y, double max_y) { + this->discrete = true; + this->min_ = min_x; + this->max_ = max_x; + this->heat_map.resize(number_of_y_pixels); + double step_x = (max_x - min_x) / (number_of_x_pixels - 1); + double step_y = (max_y - min_y) / (number_of_y_pixels - 1); + + int num_pts = interval.size(); + for (size_t i = 0; i < number_of_y_pixels; i++) { + double y = min_y + i * step_y; + this->heat_map[i].reserve(number_of_x_pixels); + for (size_t j = 0; j < number_of_x_pixels; j++) { + double x = min_x + j * step_x; + std::pair<double, double> grid_point(x, y); + double pixel_value = 0; + for (int k = 0; k < num_pts; k++) pixel_value += this->f(interval[k]) * kernel(interval[k], grid_point); + this->heat_map[i].push_back(pixel_value); + } + } + this->set_up_parameters_for_basic_classes(); +} + +template <typename Scalling_of_kernels> +Persistence_heat_maps<Scalling_of_kernels>::Persistence_heat_maps( + const std::vector<std::pair<double, double> >& interval, + const std::function<double(std::pair<double, double>, std::pair<double, double>)>& kernel) { + this->discrete = false; + this->interval = interval; + this->kernel = kernel; + int num_pts = this->interval.size(); + for (int i = 0; i < num_pts; i++) this->weights.push_back(this->f(this->interval[i])); + this->set_up_parameters_for_basic_classes(); +} + // if min_ == max_, then the program is requested to set up the values itself based on persistence intervals template <typename Scalling_of_kernels> void Persistence_heat_maps<Scalling_of_kernels>::construct(const std::vector<std::pair<double, double> >& intervals_, @@ -826,13 +888,18 @@ void Persistence_heat_maps<Scalling_of_kernels>::load_from_file(const char* file // Concretizations of virtual methods: template <typename Scalling_of_kernels> std::vector<double> Persistence_heat_maps<Scalling_of_kernels>::vectorize(int number_of_function) const { + std::vector<double> result; + if (!discrete) { + std::cout << "No vectorize method in case of infinite dimensional vectorization" << std::endl; + return result; + } + // convert this->heat_map into one large vector: size_t size_of_result = 0; for (size_t i = 0; i != this->heat_map.size(); ++i) { size_of_result += this->heat_map[i].size(); } - std::vector<double> result; result.reserve(size_of_result); for (size_t i = 0; i != this->heat_map.size(); ++i) { @@ -846,34 +913,39 @@ std::vector<double> Persistence_heat_maps<Scalling_of_kernels>::vectorize(int nu template <typename Scalling_of_kernels> double Persistence_heat_maps<Scalling_of_kernels>::distance(const Persistence_heat_maps& second, double power) const { - // first we need to check if (*this) and second are defined on the same domain and have the same dimensions: - if (!this->check_if_the_same(second)) { - std::cerr << "The persistence images are of non compatible sizes. We cannot therefore compute distance between " - "them. The program will now terminate"; - throw "The persistence images are of non compatible sizes. The program will now terminate"; - } + if (this->discrete) { + // first we need to check if (*this) and second are defined on the same domain and have the same dimensions: + if (!this->check_if_the_same(second)) { + std::cerr << "The persistence images are of non compatible sizes. We cannot therefore compute distance between " + "them. The program will now terminate"; + throw "The persistence images are of non compatible sizes. The program will now terminate"; + } - // if we are here, we know that the two persistence images are defined on the same domain, so we can start computing - // their distances: + // if we are here, we know that the two persistence images are defined on the same domain, so we can start + // computing their distances: - double distance = 0; - if (power < std::numeric_limits<double>::max()) { - for (size_t i = 0; i != this->heat_map.size(); ++i) { - for (size_t j = 0; j != this->heat_map[i].size(); ++j) { - distance += pow(fabs(this->heat_map[i][j] - second.heat_map[i][j]), power); + double distance = 0; + if (power < std::numeric_limits<double>::max()) { + for (size_t i = 0; i != this->heat_map.size(); ++i) { + for (size_t j = 0; j != this->heat_map[i].size(); ++j) { + distance += std::pow(std::fabs(this->heat_map[i][j] - second.heat_map[i][j]), power); + } } - } - } else { - // in this case, we compute max norm distance - for (size_t i = 0; i != this->heat_map.size(); ++i) { - for (size_t j = 0; j != this->heat_map[i].size(); ++j) { - if (distance < fabs(this->heat_map[i][j] - second.heat_map[i][j])) { - distance = fabs(this->heat_map[i][j] - second.heat_map[i][j]); + } else { + // in this case, we compute max norm distance + for (size_t i = 0; i != this->heat_map.size(); ++i) { + for (size_t j = 0; j != this->heat_map[i].size(); ++j) { + if (distance < std::fabs(this->heat_map[i][j] - second.heat_map[i][j])) { + distance = std::fabs(this->heat_map[i][j] - second.heat_map[i][j]); + } } } } + return distance; + } else { + return std::sqrt(this->compute_scalar_product(*this) + second.compute_scalar_product(second) - + 2 * this->compute_scalar_product(second)); } - return distance; } template <typename Scalling_of_kernels> @@ -895,22 +967,36 @@ void Persistence_heat_maps<Scalling_of_kernels>::compute_average( template <typename Scalling_of_kernels> double Persistence_heat_maps<Scalling_of_kernels>::compute_scalar_product(const Persistence_heat_maps& second) const { - // first we need to check if (*this) and second are defined on the same domain and have the same dimensions: - if (!this->check_if_the_same(second)) { - std::cerr << "The persistence images are of non compatible sizes. We cannot therefore compute distance between " - "them. The program will now terminate"; - throw "The persistence images are of non compatible sizes. The program will now terminate"; - } + if (discrete) { + // first we need to check if (*this) and second are defined on the same domain and have the same dimensions: + if (!this->check_if_the_same(second)) { + std::cerr << "The persistence images are of non compatible sizes. We cannot therefore compute distance between " + "them. The program will now terminate"; + throw "The persistence images are of non compatible sizes. The program will now terminate"; + } - // if we are here, we know that the two persistence images are defined on the same domain, so we can start computing - // their scalar product: - double scalar_prod = 0; - for (size_t i = 0; i != this->heat_map.size(); ++i) { - for (size_t j = 0; j != this->heat_map[i].size(); ++j) { - scalar_prod += this->heat_map[i][j] * second.heat_map[i][j]; + // if we are here, we know that the two persistence images are defined on the same domain, so we can start computing + // their scalar product: + double scalar_prod = 0; + for (size_t i = 0; i != this->heat_map.size(); ++i) { + for (size_t j = 0; j != this->heat_map[i].size(); ++j) { + scalar_prod += this->heat_map[i][j] * second.heat_map[i][j]; + } + } + return scalar_prod; + } else { + int num_pts1 = this->interval.size(); + int num_pts2 = second.interval.size(); + double kernel_val = 0; + for (int i = 0; i < num_pts1; i++) { + std::pair<double, double> pi = this->interval[i]; + for (int j = 0; j < num_pts2; j++) { + std::pair<double, double> pj = second.interval[j]; + kernel_val += this->weights[i] * second.weights[j] * this->kernel(pi, pj); + } } + return kernel_val; } - return scalar_prod; } } // namespace Persistence_representations diff --git a/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h new file mode 100644 index 00000000..18165c5f --- /dev/null +++ b/src/Persistence_representations/include/gudhi/Sliced_Wasserstein.h @@ -0,0 +1,392 @@ +/* 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): Mathieu Carriere + * + * 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 SLICED_WASSERSTEIN_H_ +#define SLICED_WASSERSTEIN_H_ + +// gudhi include +#include <gudhi/read_persistence_from_file.h> +#include <gudhi/common_persistence_representations.h> +#include <gudhi/Debug_utils.h> + +#include <vector> // for std::vector<> +#include <utility> // for std::pair<>, std::move +#include <algorithm> // for std::sort, std::max, std::merge +#include <cmath> // for std::abs, std::sqrt +#include <stdexcept> // for std::invalid_argument +#include <random> // for std::random_device + +namespace Gudhi { +namespace Persistence_representations { + +/** + * \class Sliced_Wasserstein gudhi/Sliced_Wasserstein.h + * \brief A class implementing the Sliced Wasserstein kernel. + * + * \ingroup Persistence_representations + * + * \details + * In this class, we compute infinite-dimensional representations of persistence diagrams by using the + * Sliced Wasserstein kernel (see \ref sec_persistence_kernels for more details on kernels). We recall that + * infinite-dimensional representations are defined implicitly, so only scalar products and distances are available for + * the representations defined in this class. + * The Sliced Wasserstein kernel is defined as a Gaussian-like kernel between persistence diagrams, where the distance + * used for comparison is the Sliced Wasserstein distance \f$SW\f$ between persistence diagrams, defined as the + * integral of the 1-norm between the sorted projections of the diagrams onto all lines passing through the origin: + * + * \f$ SW(D_1,D_2)=\int_{\theta\in\mathbb{S}}\,\|\pi_\theta(D_1\cup\pi_\Delta(D_2))-\pi_\theta(D_2\cup\pi_\Delta(D_1))\ + * |_1{\rm d}\theta\f$, + * + * where \f$\pi_\theta\f$ is the projection onto the line defined with angle \f$\theta\f$ in the unit circle + * \f$\mathbb{S}\f$, and \f$\pi_\Delta\f$ is the projection onto the diagonal. + * Assuming that the diagrams are in general position (i.e. there is no collinear triple), the integral can be computed + * exactly in \f$O(n^2{\rm log}(n))\f$ time, where \f$n\f$ is the number of points in the diagrams. We provide two + * approximations of the integral: one in which we slightly perturb the diagram points so that they are in general + * position, and another in which we approximate the integral by sampling \f$N\f$ lines in the circle in + * \f$O(Nn{\rm log}(n))\f$ time. The Sliced Wasserstein Kernel is then computed as: + * + * \f$ k(D_1,D_2) = {\rm exp}\left(-\frac{SW(D_1,D_2)}{2\sigma^2}\right).\f$ + * + * The first method is usually much more accurate but also + * much slower. For more details, please see \cite pmlr-v70-carriere17a . + * + **/ + +class Sliced_Wasserstein { + protected: + Persistence_diagram diagram; + int approx; + double sigma; + std::vector<std::vector<double> > projections, projections_diagonal; + + // ********************************** + // Utils. + // ********************************** + + void build_rep() { + if (approx > 0) { + double step = pi / this->approx; + int n = diagram.size(); + + for (int i = 0; i < this->approx; i++) { + std::vector<double> l, l_diag; + for (int j = 0; j < n; j++) { + double px = diagram[j].first; + double py = diagram[j].second; + double proj_diag = (px + py) / 2; + + l.push_back(px * cos(-pi / 2 + i * step) + py * sin(-pi / 2 + i * step)); + l_diag.push_back(proj_diag * cos(-pi / 2 + i * step) + proj_diag * sin(-pi / 2 + i * step)); + } + + std::sort(l.begin(), l.end()); + std::sort(l_diag.begin(), l_diag.end()); + projections.push_back(std::move(l)); + projections_diagonal.push_back(std::move(l_diag)); + } + + diagram.clear(); + } + } + + // Compute the angle formed by two points of a PD + double compute_angle(const Persistence_diagram& diag, int i, int j) const { + if (diag[i].second == diag[j].second) + return pi / 2; + else + return atan((diag[j].first - diag[i].first) / (diag[i].second - diag[j].second)); + } + + // Compute the integral of |cos()| between alpha and beta, valid only if alpha is in [-pi,pi] and beta-alpha is in + // [0,pi] + double compute_int_cos(double alpha, double beta) const { + double res = 0; + if (alpha >= 0 && alpha <= pi) { + if (cos(alpha) >= 0) { + if (pi / 2 <= beta) { + res = 2 - sin(alpha) - sin(beta); + } else { + res = sin(beta) - sin(alpha); + } + } else { + if (1.5 * pi <= beta) { + res = 2 + sin(alpha) + sin(beta); + } else { + res = sin(alpha) - sin(beta); + } + } + } + if (alpha >= -pi && alpha <= 0) { + if (cos(alpha) <= 0) { + if (-pi / 2 <= beta) { + res = 2 + sin(alpha) + sin(beta); + } else { + res = sin(alpha) - sin(beta); + } + } else { + if (pi / 2 <= beta) { + res = 2 - sin(alpha) - sin(beta); + } else { + res = sin(beta) - sin(alpha); + } + } + } + return res; + } + + double compute_int(double theta1, double theta2, int p, int q, const Persistence_diagram& diag1, + const Persistence_diagram& diag2) const { + double norm = std::sqrt((diag1[p].first - diag2[q].first) * (diag1[p].first - diag2[q].first) + + (diag1[p].second - diag2[q].second) * (diag1[p].second - diag2[q].second)); + double angle1; + if (diag1[p].first == diag2[q].first) + angle1 = theta1 - pi / 2; + else + angle1 = theta1 - atan((diag1[p].second - diag2[q].second) / (diag1[p].first - diag2[q].first)); + double angle2 = angle1 + theta2 - theta1; + double integral = compute_int_cos(angle1, angle2); + return norm * integral; + } + + // Evaluation of the Sliced Wasserstein Distance between a pair of diagrams. + double compute_sliced_wasserstein_distance(const Sliced_Wasserstein& second) const { + GUDHI_CHECK(this->approx == second.approx, + std::invalid_argument("Error: different approx values for representations")); + + Persistence_diagram diagram1 = this->diagram; + Persistence_diagram diagram2 = second.diagram; + double sw = 0; + + if (this->approx == -1) { + // Add projections onto diagonal. + int n1, n2; + n1 = diagram1.size(); + n2 = diagram2.size(); + double min_ordinate = std::numeric_limits<double>::max(); + double min_abscissa = std::numeric_limits<double>::max(); + double max_ordinate = std::numeric_limits<double>::lowest(); + double max_abscissa = std::numeric_limits<double>::lowest(); + for (int i = 0; i < n2; i++) { + min_ordinate = std::min(min_ordinate, diagram2[i].second); + min_abscissa = std::min(min_abscissa, diagram2[i].first); + max_ordinate = std::max(max_ordinate, diagram2[i].second); + max_abscissa = std::max(max_abscissa, diagram2[i].first); + diagram1.emplace_back((diagram2[i].first + diagram2[i].second) / 2, + (diagram2[i].first + diagram2[i].second) / 2); + } + for (int i = 0; i < n1; i++) { + min_ordinate = std::min(min_ordinate, diagram1[i].second); + min_abscissa = std::min(min_abscissa, diagram1[i].first); + max_ordinate = std::max(max_ordinate, diagram1[i].second); + max_abscissa = std::max(max_abscissa, diagram1[i].first); + diagram2.emplace_back((diagram1[i].first + diagram1[i].second) / 2, + (diagram1[i].first + diagram1[i].second) / 2); + } + int num_pts_dgm = diagram1.size(); + + // Slightly perturb the points so that the PDs are in generic positions. + double epsilon = 0.0001; + double thresh_y = (max_ordinate - min_ordinate) * epsilon; + double thresh_x = (max_abscissa - min_abscissa) * epsilon; + std::random_device rd; + std::default_random_engine re(rd()); + std::uniform_real_distribution<double> uni(-1, 1); + for (int i = 0; i < num_pts_dgm; i++) { + double u = uni(re); + diagram1[i].first += u * thresh_x; + diagram1[i].second += u * thresh_y; + diagram2[i].first += u * thresh_x; + diagram2[i].second += u * thresh_y; + } + + // Compute all angles in both PDs. + std::vector<std::pair<double, std::pair<int, int> > > angles1, angles2; + for (int i = 0; i < num_pts_dgm; i++) { + for (int j = i + 1; j < num_pts_dgm; j++) { + double theta1 = compute_angle(diagram1, i, j); + double theta2 = compute_angle(diagram2, i, j); + angles1.emplace_back(theta1, std::pair<int, int>(i, j)); + angles2.emplace_back(theta2, std::pair<int, int>(i, j)); + } + } + + // Sort angles. + std::sort(angles1.begin(), angles1.end(), + [](const std::pair<double, std::pair<int, int> >& p1, + const std::pair<double, std::pair<int, int> >& p2) { return (p1.first < p2.first); }); + std::sort(angles2.begin(), angles2.end(), + [](const std::pair<double, std::pair<int, int> >& p1, + const std::pair<double, std::pair<int, int> >& p2) { return (p1.first < p2.first); }); + + // Initialize orders of the points of both PDs (given by ordinates when theta = -pi/2). + std::vector<int> orderp1, orderp2; + for (int i = 0; i < num_pts_dgm; i++) { + orderp1.push_back(i); + orderp2.push_back(i); + } + std::sort(orderp1.begin(), orderp1.end(), [&](int i, int j) { + if (diagram1[i].second != diagram1[j].second) + return (diagram1[i].second < diagram1[j].second); + else + return (diagram1[i].first > diagram1[j].first); + }); + std::sort(orderp2.begin(), orderp2.end(), [&](int i, int j) { + if (diagram2[i].second != diagram2[j].second) + return (diagram2[i].second < diagram2[j].second); + else + return (diagram2[i].first > diagram2[j].first); + }); + + // Find the inverses of the orders. + std::vector<int> order1(num_pts_dgm); + std::vector<int> order2(num_pts_dgm); + for (int i = 0; i < num_pts_dgm; i++) { + order1[orderp1[i]] = i; + order2[orderp2[i]] = i; + } + + // Record all inversions of points in the orders as theta varies along the positive half-disk. + std::vector<std::vector<std::pair<int, double> > > anglePerm1(num_pts_dgm); + std::vector<std::vector<std::pair<int, double> > > anglePerm2(num_pts_dgm); + + int m1 = angles1.size(); + for (int i = 0; i < m1; i++) { + double theta = angles1[i].first; + int p = angles1[i].second.first; + int q = angles1[i].second.second; + anglePerm1[order1[p]].emplace_back(p, theta); + anglePerm1[order1[q]].emplace_back(q, theta); + int a = order1[p]; + int b = order1[q]; + order1[p] = b; + order1[q] = a; + } + + int m2 = angles2.size(); + for (int i = 0; i < m2; i++) { + double theta = angles2[i].first; + int p = angles2[i].second.first; + int q = angles2[i].second.second; + anglePerm2[order2[p]].emplace_back(p, theta); + anglePerm2[order2[q]].emplace_back(q, theta); + int a = order2[p]; + int b = order2[q]; + order2[p] = b; + order2[q] = a; + } + + for (int i = 0; i < num_pts_dgm; i++) { + anglePerm1[order1[i]].emplace_back(i, pi / 2); + anglePerm2[order2[i]].emplace_back(i, pi / 2); + } + + // Compute the SW distance with the list of inversions. + for (int i = 0; i < num_pts_dgm; i++) { + std::vector<std::pair<int, double> > u, v; + u = anglePerm1[i]; + v = anglePerm2[i]; + double theta1, theta2; + theta1 = -pi / 2; + unsigned int ku, kv; + ku = 0; + kv = 0; + theta2 = std::min(u[ku].second, v[kv].second); + while (theta1 != pi / 2) { + if (diagram1[u[ku].first].first != diagram2[v[kv].first].first || + diagram1[u[ku].first].second != diagram2[v[kv].first].second) + if (theta1 != theta2) sw += compute_int(theta1, theta2, u[ku].first, v[kv].first, diagram1, diagram2); + theta1 = theta2; + if ((theta2 == u[ku].second) && ku < u.size() - 1) ku++; + if ((theta2 == v[kv].second) && kv < v.size() - 1) kv++; + theta2 = std::min(u[ku].second, v[kv].second); + } + } + } else { + double step = pi / this->approx; + std::vector<double> v1, v2; + for (int i = 0; i < this->approx; i++) { + v1.clear(); + v2.clear(); + std::merge(this->projections[i].begin(), this->projections[i].end(), second.projections_diagonal[i].begin(), + second.projections_diagonal[i].end(), std::back_inserter(v1)); + std::merge(second.projections[i].begin(), second.projections[i].end(), this->projections_diagonal[i].begin(), + this->projections_diagonal[i].end(), std::back_inserter(v2)); + + int n = v1.size(); + double f = 0; + for (int j = 0; j < n; j++) f += std::abs(v1[j] - v2[j]); + sw += f * step; + } + } + + return sw / pi; + } + + public: + /** \brief Sliced Wasserstein kernel constructor. + * \implements Topological_data_with_distances, Real_valued_topological_data, Topological_data_with_scalar_product + * \ingroup Sliced_Wasserstein + * + * @param[in] _diagram persistence diagram. + * @param[in] _sigma bandwidth parameter. + * @param[in] _approx number of directions used to approximate the integral in the Sliced Wasserstein distance, set + * to -1 for random perturbation. If positive, then projections of the diagram points on all + * directions are stored in memory to reduce computation time. + * + */ + Sliced_Wasserstein(const Persistence_diagram& _diagram, double _sigma = 1.0, int _approx = 10) + : diagram(_diagram), approx(_approx), sigma(_sigma) { + build_rep(); + } + + /** \brief Evaluation of the kernel on a pair of diagrams. + * \ingroup Sliced_Wasserstein + * + * @pre approx and sigma attributes need to be the same for both instances. + * @param[in] second other instance of class Sliced_Wasserstein. + * + */ + double compute_scalar_product(const Sliced_Wasserstein& second) const { + GUDHI_CHECK(this->sigma == second.sigma, + std::invalid_argument("Error: different sigma values for representations")); + return std::exp(-compute_sliced_wasserstein_distance(second) / (2 * this->sigma * this->sigma)); + } + + /** \brief Evaluation of the distance between images of diagrams in the Hilbert space of the kernel. + * \ingroup Sliced_Wasserstein + * + * @pre approx and sigma attributes need to be the same for both instances. + * @param[in] second other instance of class Sliced_Wasserstein. + * + */ + double distance(const Sliced_Wasserstein& second) const { + GUDHI_CHECK(this->sigma == second.sigma, + std::invalid_argument("Error: different sigma values for representations")); + return std::sqrt(this->compute_scalar_product(*this) + second.compute_scalar_product(second) - + 2 * this->compute_scalar_product(second)); + } + +}; // class Sliced_Wasserstein +} // namespace Persistence_representations +} // namespace Gudhi + +#endif // SLICED_WASSERSTEIN_H_ diff --git a/src/Persistence_representations/include/gudhi/common_persistence_representations.h b/src/Persistence_representations/include/gudhi/common_persistence_representations.h index 3d03f1f6..6fed019a 100644 --- a/src/Persistence_representations/include/gudhi/common_persistence_representations.h +++ b/src/Persistence_representations/include/gudhi/common_persistence_representations.h @@ -26,17 +26,28 @@ #include <utility> #include <string> #include <cmath> +#include <boost/math/constants/constants.hpp> + + namespace Gudhi { namespace Persistence_representations { // this file contain an implementation of some common procedures used in Persistence_representations. +static constexpr double pi = boost::math::constants::pi<double>(); + + +/** + * In this module, we use the name Persistence_diagram for the representation of a diagram in a vector of pairs of two double. + */ +using Persistence_diagram = std::vector<std::pair<double, double> >; + // double epsi = std::numeric_limits<double>::epsilon(); double epsi = 0.000005; /** * A procedure used to compare doubles. Typically given two doubles A and B, comparing A == B is not good idea. In this - *case, we use the procedure almostEqual with the epsi defined at + * case, we use the procedure almostEqual with the epsi defined at * the top of the file. Setting up the epsi gives the user a tolerance on what should be consider equal. **/ inline bool almost_equal(double a, double b) { @@ -53,8 +64,7 @@ double birth_plus_deaths(std::pair<double, double> a) { return a.first + a.secon // landscapes /** - * Given two points in R^2, the procedure compute the parameters A and B of the line y = Ax + B that crosses those two - *points. + * Given two points in R^2, the procedure compute the parameters A and B of the line y = Ax + B that crosses those two points. **/ std::pair<double, double> compute_parameters_of_a_line(std::pair<double, double> p1, std::pair<double, double> p2) { double a = (p2.second - p1.second) / (p2.first - p1.first); @@ -64,8 +74,7 @@ std::pair<double, double> compute_parameters_of_a_line(std::pair<double, double> // landscapes /** - * This procedure given two points which lies on the opposite sides of x axis, compute x for which the line connecting - *those two points crosses x axis. + * This procedure given two points which lies on the opposite sides of x axis, compute x for which the line connecting those two points crosses x axis. **/ double find_zero_of_a_line_segment_between_those_two_points(std::pair<double, double> p1, std::pair<double, double> p2) { @@ -89,8 +98,7 @@ double find_zero_of_a_line_segment_between_those_two_points(std::pair<double, do // landscapes /** * This method provides a comparison of points that is used in construction of persistence landscapes. The ordering is - *lexicographical for the first coordinate, and reverse-lexicographical for the - * second coordinate. + * lexicographical for the first coordinate, and reverse-lexicographical for the second coordinate. **/ bool compare_points_sorting(std::pair<double, double> f, std::pair<double, double> s) { if (f.first < s.first) { diff --git a/src/Persistence_representations/test/CMakeLists.txt b/src/Persistence_representations/test/CMakeLists.txt index 5e2b6910..fb650485 100644 --- a/src/Persistence_representations/test/CMakeLists.txt +++ b/src/Persistence_representations/test/CMakeLists.txt @@ -34,6 +34,11 @@ target_link_libraries(Read_persistence_from_file_test_unit ${Boost_UNIT_TEST_FRA gudhi_add_coverage_test(Read_persistence_from_file_test_unit) +add_executable ( kernels_unit kernels.cpp ) +target_link_libraries(kernels_unit ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + +gudhi_add_coverage_test(kernels_unit) + if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1) add_executable (Persistence_intervals_with_distances_test_unit persistence_intervals_with_distances_test.cpp ) target_link_libraries(Persistence_intervals_with_distances_test_unit ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) diff --git a/src/Persistence_representations/test/kernels.cpp b/src/Persistence_representations/test/kernels.cpp new file mode 100644 index 00000000..b8d02d4c --- /dev/null +++ b/src/Persistence_representations/test/kernels.cpp @@ -0,0 +1,61 @@ +/* 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): Mathieu Carrière + * + * 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/>. + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "kernel" + +#include <boost/test/unit_test.hpp> +#include <cmath> // float comparison +#include <limits> +#include <string> +#include <vector> +#include <algorithm> // std::max +#include <gudhi/Persistence_heat_maps.h> +#include <gudhi/common_persistence_representations.h> +#include <gudhi/Sliced_Wasserstein.h> +#include <gudhi/distance_functions.h> +#include <gudhi/reader_utils.h> + +using constant_scaling_function = Gudhi::Persistence_representations::constant_scaling_function; +using SW = Gudhi::Persistence_representations::Sliced_Wasserstein; +using PWG = Gudhi::Persistence_representations::Persistence_heat_maps<constant_scaling_function>; +using Persistence_diagram = std::vector<std::pair<double,double> >; + +std::function<double(std::pair<double, double>, std::pair<double, double>)> Gaussian_function(double sigma){ + return [=](std::pair<double, double> p, std::pair<double, double> q){ + return (1/std::sqrt(2*Gudhi::Persistence_representations::pi)*sigma) * std::exp( -( (p.first-q.first)*(p.first-q.first) + (p.second-q.second)*(p.second-q.second) )/(2*sigma) ); + }; +} + +BOOST_AUTO_TEST_CASE(check_PWG) { + Persistence_diagram v1, v2; v1.emplace_back(0,1); v2.emplace_back(0,2); + PWG pwg1(v1, Gaussian_function(1.0)); + PWG pwg2(v2, Gaussian_function(1.0)); + BOOST_CHECK(std::abs(pwg1.compute_scalar_product(pwg2) - std::exp(-0.5)/(std::sqrt(2*Gudhi::Persistence_representations::pi))) <= 1e-3); +} + +BOOST_AUTO_TEST_CASE(check_SW) { + Persistence_diagram v1, v2; v1.emplace_back(0,1); v2.emplace_back(0,2); + SW sw1(v1, 1.0, 100); SW swex1(v1, 1.0, -1); + SW sw2(v2, 1.0, 100); SW swex2(v2, 1.0, -1); + BOOST_CHECK(std::abs(sw1.compute_scalar_product(sw2) - swex1.compute_scalar_product(swex2)) <= 1e-1); +} diff --git a/src/Persistent_cohomology/concept/FilteredComplex.h b/src/Persistent_cohomology/concept/FilteredComplex.h index 62b9002f..7eb01b01 100644 --- a/src/Persistent_cohomology/concept/FilteredComplex.h +++ b/src/Persistent_cohomology/concept/FilteredComplex.h @@ -27,7 +27,7 @@ */ struct FilteredComplex { -/** Handle to specify a simplex. */ +/** \brief Handle to specify a simplex. */ typedef unspecified Simplex_handle; /** \brief Type for the value of the filtration function. * @@ -39,8 +39,7 @@ struct FilteredComplex * is model of IndexingTag. */ typedef unspecified Indexing_tag; -/** Returns a Simplex_handle that is different from all simplex handles - * of the simplices. */ +/** \brief Returns a Simplex_handle that is different from all simplex handles of the simplices. */ Simplex_handle null_simplex(); /** \brief Returns the number of simplices in the complex. * @@ -58,22 +57,19 @@ struct FilteredComplex * * This is only called on valid indices. */ Simplex_handle simplex ( size_t idx ); -/** \brief Iterator on the simplices belonging to the - * boundary of a simplex. +/** \brief Iterator on the simplices belonging to the boundary of a simplex. * * <CODE>value_type</CODE> must be 'Simplex_handle'. */ typedef unspecified Boundary_simplex_iterator; -/** \brief Range giving access to the simplices in the boundary of - * a simplex. +/** \brief Range giving access to the simplices in the boundary of a simplex. * * .begin() and .end() return type Boundary_simplex_iterator. */ typedef unspecified Boundary_simplex_range; -/** \brief Returns a range giving access to all simplices of the - * boundary of a simplex, i.e. - * the set of codimension 1 subsimplices of the Simplex. +/** \brief Returns a range giving access to all simplices of the boundary of a simplex, i.e. the set of codimension 1 + * subsimplices of the Simplex. * * If the simplex is \f$[v_0, \cdots ,v_d]\f$, with canonical orientation * induced by \f$ v_0 < \cdots < v_d \f$, the iterator enumerates the @@ -84,19 +80,16 @@ typedef unspecified Boundary_simplex_range; * gives the chains corresponding to the boundary of the simplex.*/ Boundary_simplex_range boundary_simplex_range(Simplex_handle sh); -/** \brief Iterator over all simplices of the complex - * in the order of the indexing scheme. +/** \brief Iterator over all simplices of the complex in the order of the indexing scheme. * * 'value_type' must be 'Simplex_handle'. */ typedef unspecified Filtration_simplex_iterator; -/** \brief Range over the simplices of the complex - * in the order of the filtration. +/** \brief Range over the simplices of the complex in the order of the filtration. * * .begin() and .end() return type Filtration_simplex_iterator.*/ typedef unspecified Filtration_simplex_range; -/** \brief Returns a range over the simplices of the complex - * in the order of the filtration. +/** \brief Returns a range over the simplices of the complex in the order of the filtration. * * .begin() and .end() return type Filtration_simplex_iterator.*/ Filtration_simplex_range filtration_simplex_range(); diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h index c51e47a5..c57174cb 100644 --- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h @@ -63,12 +63,21 @@ namespace persistent_cohomology { template<class FilteredComplex, class CoefficientField> class Persistent_cohomology { public: - typedef FilteredComplex Complex_ds; // Data attached to each simplex to interface with a Property Map. - typedef typename Complex_ds::Simplex_key Simplex_key; - typedef typename Complex_ds::Simplex_handle Simplex_handle; - typedef typename Complex_ds::Filtration_value Filtration_value; + + /** \brief Data stored for each simplex. */ + typedef typename FilteredComplex::Simplex_key Simplex_key; + /** \brief Handle to specify a simplex. */ + typedef typename FilteredComplex::Simplex_handle Simplex_handle; + /** \brief Type for the value of the filtration function. */ + typedef typename FilteredComplex::Filtration_value Filtration_value; + /** \brief Type of element of the field. */ typedef typename CoefficientField::Element Arith_element; + /** \brief Type for birth and death FilteredComplex::Simplex_handle. + * The Arith_element field is used for the multi-field framework. */ + typedef std::tuple<Simplex_handle, Simplex_handle, Arith_element> Persistent_interval; + + private: // Compressed Annotation Matrix types: // Column type typedef Persistent_cohomology_column<Simplex_key, Arith_element> Column; // contains 1 set_hook @@ -83,16 +92,15 @@ class Persistent_cohomology { boost::intrusive::constant_time_size<false> > Cam; // Sparse column type for the annotation of the boundary of an element. typedef std::vector<std::pair<Simplex_key, Arith_element> > A_ds_type; - // Persistent interval type. The Arith_element field is used for the multi-field framework. - typedef std::tuple<Simplex_handle, Simplex_handle, Arith_element> Persistent_interval; + public: /** \brief Initializes the Persistent_cohomology class. * * @param[in] cpx Complex for which the persistent homology is computed. * cpx is a model of FilteredComplex * @exception std::out_of_range In case the number of simplices is more than Simplex_key type numeric limit. */ - explicit Persistent_cohomology(Complex_ds& cpx) + explicit Persistent_cohomology(FilteredComplex& cpx) : cpx_(&cpx), dim_max_(cpx.dimension()), // upper bound on the dimension of the simplices coeff_field_(), // initialize the field coefficient structure. @@ -128,7 +136,7 @@ class Persistent_cohomology { * @param[in] persistence_dim_max if true, the persistent homology for the maximal dimension in the * complex is computed. If false, it is ignored. Default is false. */ - Persistent_cohomology(Complex_ds& cpx, bool persistence_dim_max) + Persistent_cohomology(FilteredComplex& cpx, bool persistence_dim_max) : Persistent_cohomology(cpx) { if (persistence_dim_max) { ++dim_max_; @@ -146,7 +154,7 @@ class Persistent_cohomology { private: struct length_interval { - length_interval(Complex_ds * cpx, Filtration_value min_length) + length_interval(FilteredComplex * cpx, Filtration_value min_length) : cpx_(cpx), min_length_(min_length) { } @@ -159,7 +167,7 @@ class Persistent_cohomology { min_length_ = new_length; } - Complex_ds * cpx_; + FilteredComplex * cpx_; Filtration_value min_length_; }; @@ -552,14 +560,14 @@ class Persistent_cohomology { * Compare two intervals by length. */ struct cmp_intervals_by_length { - explicit cmp_intervals_by_length(Complex_ds * sc) + explicit cmp_intervals_by_length(FilteredComplex * sc) : sc_(sc) { } bool operator()(const Persistent_interval & p1, const Persistent_interval & p2) { return (sc_->filtration(get < 1 > (p1)) - sc_->filtration(get < 0 > (p1)) > sc_->filtration(get < 1 > (p2)) - sc_->filtration(get < 0 > (p2))); } - Complex_ds * sc_; + FilteredComplex * sc_; }; public: @@ -690,9 +698,8 @@ class Persistent_cohomology { return betti_number; } - /** @brief Returns the persistent pairs. - * @return Persistent pairs - * + /** @brief Returns a list of persistence birth and death FilteredComplex::Simplex_handle pairs. + * @return A list of Persistent_cohomology::Persistent_interval */ const std::vector<Persistent_interval>& get_persistent_pairs() const { return persistent_pairs_; @@ -733,7 +740,7 @@ class Persistent_cohomology { }; public: - Complex_ds * cpx_; + FilteredComplex * cpx_; int dim_max_; CoefficientField coeff_field_; size_t num_simplices_; diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 4b18651c..343ed472 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -1330,7 +1330,7 @@ class Simplex_tree { public: /** \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. + * @return True if any filtration value was modified, false if the filtration was already non-decreasing. * \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 * complex has changed , please call `initialize_filtration()` to recompute it. diff --git a/src/common/doc/MathJax.js b/src/common/doc/MathJax.js deleted file mode 100644 index 35e1994e..00000000 --- a/src/common/doc/MathJax.js +++ /dev/null @@ -1,53 +0,0 @@ -(function () { - var newMathJax = 'https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js'; - var oldMathJax = 'cdn.mathjax.org/mathjax/latest/MathJax.js'; - - var replaceScript = function (script, src) { - // - // Make redirected script - // - var newScript = document.createElement('script'); - newScript.src = newMathJax + src.replace(/.*?(\?|$)/, '$1'); - // - // Move onload and onerror handlers to new script - // - newScript.onload = script.onload; - newScript.onerror = script.onerror; - script.onload = script.onerror = null; - // - // Move any content (old-style configuration scripts) - // - while (script.firstChild) newScript.appendChild(script.firstChild); - // - // Copy script id - // - if (script.id != null) newScript.id = script.id; - // - // Replace original script with new one - // - script.parentNode.replaceChild(newScript, script); - // - // Issue a console warning - // - console.warn('WARNING: cdn.mathjax.org has been retired. Check https://www.mathjax.org/cdn-shutting-down/ for migration tips.') - } - - if (document.currentScript) { - var script = document.currentScript; - replaceScript(script, script.src); - } else { - // - // Look for current script by searching for one with the right source - // - var n = oldMathJax.length; - var scripts = document.getElementsByTagName('script'); - for (var i = 0; i < scripts.length; i++) { - var script = scripts[i]; - var src = (script.src || '').replace(/.*?:\/\//,''); - if (src.substr(0, n) === oldMathJax) { - replaceScript(script, src); - break; - } - } - } -})();
\ No newline at end of file diff --git a/src/common/doc/file_formats.h b/src/common/doc/file_formats.h index 23214e25..235296d3 100644 --- a/src/common/doc/file_formats.h +++ b/src/common/doc/file_formats.h @@ -29,6 +29,32 @@ namespace Gudhi { \tableofcontents + \section FileFormatsOFF OFF file format + + OFF files must be conform to format described here: http://www.geomview.org/docs/html/OFF.html + + OFF files are mainly used as point cloud inputs. Here is an example of 7 points in a 3-dimensional space. As edges and + faces are not used for point set, there is no need to specify them (just set their numbers to 0): + + \include points/alphacomplexdoc.off + + For dimensions bigger than 3, the dimension can be set like here: + \verbatim + # Dimension is no more 3 + nOFF + # dimension 4 - 7 vertices - 0 face - 0 edge + 4 7 0 0 + # Point set: + 1.0 1.0 0.0 0.0 + 7.0 0.0 0.0 0.0 + 4.0 6.0 0.0 0.0 + 9.0 6.0 0.0 0.0 + 0.0 14.0 0.0 0.0 + 2.0 19.0 0.0 0.0 + 9.0 17.0 0.0 0.0 + \endverbatim + + \section FileFormatsPers Persistence Diagram Such a file, whose extension is usually `.pers`, contains a list of persistence intervals.<br> diff --git a/src/common/doc/installation.h b/src/common/doc/installation.h index 8fb8b330..5d581b08 100644 --- a/src/common/doc/installation.h +++ b/src/common/doc/installation.h @@ -44,7 +44,7 @@ make doxygen \endverbatim * * \subsection helloworld Hello world ! - * The <a target="_blank" href="https://gitlab.inria.fr/GUDHI/hello-gudhi-world">Hello world for GUDHI</a> + * The <a target="_blank" href="https://github.com/GUDHI/hello-gudhi-world">Hello world for GUDHI</a> * project is an example to help developers to make their own C++ project on top of the GUDHI library. * * \section optionallibrary Optional third-party library diff --git a/src/common/example/vectordoubleoffreader_result.txt b/src/common/example/vectordoubleoffreader_result.txt index 1deb8dbd..b399425a 100644 --- a/src/common/example/vectordoubleoffreader_result.txt +++ b/src/common/example/vectordoubleoffreader_result.txt @@ -1,7 +1,7 @@ -Point[0] = 1 1 -Point[1] = 7 0 -Point[2] = 4 6 -Point[3] = 9 6 -Point[4] = 0 14 -Point[5] = 2 19 -Point[6] = 9 17 +Point[0] = 1 1 0 +Point[1] = 7 0 0 +Point[2] = 4 6 0 +Point[3] = 9 6 0 +Point[4] = 0 14 0 +Point[5] = 2 19 0 +Point[6] = 9 17 0 diff --git a/src/common/include/gudhi/Off_reader.h b/src/common/include/gudhi/Off_reader.h index 05a1e145..fc951fe7 100644 --- a/src/common/include/gudhi/Off_reader.h +++ b/src/common/include/gudhi/Off_reader.h @@ -37,8 +37,7 @@ namespace Gudhi { /** \brief OFF file reader top class visitor. * - * OFF file must be conform to format described here : - * http://www.geomview.org/docs/html/OFF.html + * OFF file must be conform to \ref FileFormatsOFF */ class Off_reader { public: diff --git a/src/common/test/test_points_off_reader.cpp b/src/common/test/test_points_off_reader.cpp index ba3bab71..e4b76ed7 100644 --- a/src/common/test/test_points_off_reader.cpp +++ b/src/common/test/test_points_off_reader.cpp @@ -44,19 +44,19 @@ BOOST_AUTO_TEST_CASE( points_doc_test ) BOOST_CHECK(point_cloud.size() == 7); std::vector<Point_d> expected_points; - std::vector<double> point = {1.0, 1.0}; + std::vector<double> point = {1.0, 1.0, 0.0}; expected_points.push_back(Point_d(point.begin(), point.end())); - point = {7.0, 0.0}; + point = {7.0, 0.0, 0.0}; expected_points.push_back(Point_d(point.begin(), point.end())); - point = {4.0, 6.0}; + point = {4.0, 6.0, 0.0}; expected_points.push_back(Point_d(point.begin(), point.end())); - point = {9.0, 6.0}; + point = {9.0, 6.0, 0.0}; expected_points.push_back(Point_d(point.begin(), point.end())); - point = {0.0, 14.0}; + point = {0.0, 14.0, 0.0}; expected_points.push_back(Point_d(point.begin(), point.end())); - point = {2.0, 19.0}; + point = {2.0, 19.0, 0.0}; expected_points.push_back(Point_d(point.begin(), point.end())); - point = {9.0, 17.0}; + point = {9.0, 17.0, 0.0}; expected_points.push_back(Point_d(point.begin(), point.end())); BOOST_CHECK(point_cloud == expected_points); diff --git a/src/cython/CMakeLists.txt b/src/cython/CMakeLists.txt index 480332d7..d4ace20e 100644 --- a/src/cython/CMakeLists.txt +++ b/src/cython/CMakeLists.txt @@ -138,6 +138,7 @@ if(PYTHONINTERP_FOUND) else() add_gudhi_cython_lib("${Boost_THREAD_LIBRARY_RELEASE}") endif() + message("** Add Boost ${Boost_LIBRARY_DIRS}") set(GUDHI_CYTHON_LIBRARY_DIRS "${GUDHI_CYTHON_LIBRARY_DIRS}'${Boost_LIBRARY_DIRS}', ") endif() # Add CGAL compilation args @@ -148,6 +149,7 @@ if(PYTHONINTERP_FOUND) add_gudhi_debug_info("CGAL version ${CGAL_VERSION}") add_gudhi_cython_lib("${CGAL_LIBRARY}") set(GUDHI_CYTHON_LIBRARY_DIRS "${GUDHI_CYTHON_LIBRARY_DIRS}'${CGAL_LIBRARIES_DIR}', ") + message("** Add CGAL ${CGAL_LIBRARIES_DIR}") # If CGAL is not header only, CGAL library may link with boost system, if(CMAKE_BUILD_TYPE MATCHES Debug) add_gudhi_cython_lib("${Boost_SYSTEM_LIBRARY_DEBUG}") @@ -155,6 +157,7 @@ if(PYTHONINTERP_FOUND) add_gudhi_cython_lib("${Boost_SYSTEM_LIBRARY_RELEASE}") endif() set(GUDHI_CYTHON_LIBRARY_DIRS "${GUDHI_CYTHON_LIBRARY_DIRS}'${Boost_LIBRARY_DIRS}', ") + message("** Add Boost ${Boost_LIBRARY_DIRS}") endif(CGAL_HEADER_ONLY) # GMP and GMPXX are not required, but if present, CGAL will link with them. if(GMP_FOUND) @@ -162,11 +165,13 @@ if(PYTHONINTERP_FOUND) set(GUDHI_CYTHON_EXTRA_COMPILE_ARGS "${GUDHI_CYTHON_EXTRA_COMPILE_ARGS}'-DCGAL_USE_GMP', ") add_gudhi_cython_lib("${GMP_LIBRARIES}") set(GUDHI_CYTHON_LIBRARY_DIRS "${GUDHI_CYTHON_LIBRARY_DIRS}'${GMP_LIBRARIES_DIR}', ") + message("** Add gmp ${GMP_LIBRARIES_DIR}") if(GMPXX_FOUND) add_gudhi_debug_info("GMPXX_LIBRARIES = ${GMPXX_LIBRARIES}") set(GUDHI_CYTHON_EXTRA_COMPILE_ARGS "${GUDHI_CYTHON_EXTRA_COMPILE_ARGS}'-DCGAL_USE_GMPXX', ") add_gudhi_cython_lib("${GMPXX_LIBRARIES}") set(GUDHI_CYTHON_LIBRARY_DIRS "${GUDHI_CYTHON_LIBRARY_DIRS}'${GMPXX_LIBRARIES_DIR}', ") + message("** Add gmpxx ${GMPXX_LIBRARIES_DIR}") endif(GMPXX_FOUND) endif(GMP_FOUND) endif(CGAL_FOUND) @@ -195,6 +200,7 @@ if(PYTHONINTERP_FOUND) add_gudhi_cython_lib("${TBB_MALLOC_RELEASE_LIBRARY}") endif() set(GUDHI_CYTHON_LIBRARY_DIRS "${GUDHI_CYTHON_LIBRARY_DIRS}'${TBB_LIBRARY_DIRS}', ") + message("** Add tbb ${TBB_LIBRARY_DIRS}") set(GUDHI_CYTHON_INCLUDE_DIRS "${GUDHI_CYTHON_INCLUDE_DIRS}'${TBB_INCLUDE_DIRS}', ") endif() diff --git a/src/cython/cython/simplex_tree.pyx b/src/cython/cython/simplex_tree.pyx index 0ab97f80..ea99c940 100644 --- a/src/cython/cython/simplex_tree.pyx +++ b/src/cython/cython/simplex_tree.pyx @@ -405,7 +405,8 @@ cdef class SimplexTree: """This function ensures that each simplex has a higher filtration value than its faces by increasing the filtration values. - :returns: The filtration modification information. + :returns: True if any filtration value was modified, + False if the filtration was already non-decreasing. :rtype: bool @@ -513,10 +514,9 @@ cdef class SimplexTree: return intervals_result def persistence_pairs(self): - """This function returns the persistence pairs of the simplicial - complex. + """This function returns a list of persistence birth and death simplices pairs. - :returns: The persistence intervals. + :returns: A list of persistence simplices intervals. :rtype: list of pair of list of int :note: persistence_pairs function requires diff --git a/src/cython/doc/fileformats.rst b/src/cython/doc/fileformats.rst index e205cc8b..345dfdba 100644 --- a/src/cython/doc/fileformats.rst +++ b/src/cython/doc/fileformats.rst @@ -5,6 +5,35 @@ File formats ############ +OFF file format +*************** + +OFF files must be conform to format described here: +http://www.geomview.org/docs/html/OFF.html + +OFF files are mainly used as point cloud inputs. Here is an example of 7 points +in a 3-dimensional space. As edges and faces are not used for point set, there +is no need to specify them (just set their numbers to 0): + +.. literalinclude:: ../../data/points/alphacomplexdoc.off + +.. centered:: ../../points/alphacomplexdoc.off + +For dimensions bigger than 3, the dimension can be set like here:: + + # Dimension is no more 3 + nOFF + # dimension 4 - 7 vertices - 0 face - 0 edge + 4 7 0 0 + # Point set: + 1.0 1.0 0.0 0.0 + 7.0 0.0 0.0 0.0 + 4.0 6.0 0.0 0.0 + 9.0 6.0 0.0 0.0 + 0.0 14.0 0.0 0.0 + 2.0 19.0 0.0 0.0 + 9.0 17.0 0.0 0.0 + Persistence Diagram ******************* diff --git a/src/cython/doc/nerve_gic_complex_user.rst b/src/cython/doc/nerve_gic_complex_user.rst index 94a2b246..9101f45d 100644 --- a/src/cython/doc/nerve_gic_complex_user.rst +++ b/src/cython/doc/nerve_gic_complex_user.rst @@ -13,8 +13,7 @@ Visualizations of the simplicial complexes can be done with either neato (from `graphviz <http://www.graphviz.org/>`_), `geomview <http://www.geomview.org/>`_, `KeplerMapper <https://github.com/MLWave/kepler-mapper>`_. -Input point clouds are assumed to be -`OFF files <http://www.geomview.org/docs/html/OFF.html>`_. +Input point clouds are assumed to be OFF files (cf. :doc:`fileformats`). Covers ------ diff --git a/src/cython/include/Simplex_tree_interface.h b/src/cython/include/Simplex_tree_interface.h index 3481eeff..ca98517d 100644 --- a/src/cython/include/Simplex_tree_interface.h +++ b/src/cython/include/Simplex_tree_interface.h @@ -45,7 +45,7 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { using Simplex_handle = typename Base::Simplex_handle; using Insertion_result = typename std::pair<Simplex_handle, bool>; using Simplex = std::vector<Vertex_handle>; - using Complex = std::vector<std::pair<Simplex, Filtration_value>>; + using Filtered_simplices = std::vector<std::pair<Simplex, Filtration_value>>; public: bool find_simplex(const Simplex& vh) { @@ -94,9 +94,9 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { Base::initialize_filtration(); } - Complex get_filtration() { + Filtered_simplices get_filtration() { Base::initialize_filtration(); - Complex filtrations; + Filtered_simplices filtrations; for (auto f_simplex : Base::filtration_simplex_range()) { Simplex simplex; for (auto vertex : Base::simplex_vertex_range(f_simplex)) { @@ -107,8 +107,8 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { return filtrations; } - Complex get_skeleton(int dimension) { - Complex skeletons; + Filtered_simplices get_skeleton(int dimension) { + Filtered_simplices skeletons; for (auto f_simplex : Base::skeleton_simplex_range(dimension)) { Simplex simplex; for (auto vertex : Base::simplex_vertex_range(f_simplex)) { @@ -119,29 +119,25 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { return skeletons; } - Complex get_star(const Simplex& simplex) { - Complex star; + Filtered_simplices get_star(const Simplex& simplex) { + Filtered_simplices star; for (auto f_simplex : Base::star_simplex_range(Base::find(simplex))) { Simplex simplex_star; for (auto vertex : Base::simplex_vertex_range(f_simplex)) { - std::cout << vertex << " "; simplex_star.insert(simplex_star.begin(), vertex); } - std::cout << std::endl; star.push_back(std::make_pair(simplex_star, Base::filtration(f_simplex))); } return star; } - Complex get_cofaces(const Simplex& simplex, int dimension) { - Complex cofaces; + Filtered_simplices get_cofaces(const Simplex& simplex, int dimension) { + Filtered_simplices cofaces; for (auto f_simplex : Base::cofaces_simplex_range(Base::find(simplex), dimension)) { Simplex simplex_coface; for (auto vertex : Base::simplex_vertex_range(f_simplex)) { - std::cout << vertex << " "; simplex_coface.insert(simplex_coface.begin(), vertex); } - std::cout << std::endl; cofaces.push_back(std::make_pair(simplex_coface, Base::filtration(f_simplex))); } return cofaces; |