diff options
| author | ROUVREAU Vincent <vincent.rouvreau@inria.fr> | 2019-06-12 08:07:22 +0200 |
|---|---|---|
| committer | ROUVREAU Vincent <vincent.rouvreau@inria.fr> | 2019-06-12 08:07:22 +0200 |
| commit | 4930211002c29b78d0e090aa482f6a5b9226ef51 (patch) | |
| tree | f13addfc5fb943a68ed7ebee84a5023f41ab0697 /src | |
| parent | 2927abe9a4b00bd291703340daaaed7bf20f3753 (diff) | |
| parent | f58f0bb2cb99076d0cd3a11ad39f3277213e3f5e (diff) | |
Merge branch 'master' into simplex_tree_insert_duplicated_vertices_fix_vincent
Diffstat (limited to 'src')
52 files changed, 1485 insertions, 704 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 1c293d1c..f80d4505 100644 --- a/src/Doxyfile.in +++ b/src/Doxyfile.in @@ -208,7 +208,7 @@ SEPARATE_MEMBER_PAGES = NO # uses this value to replace tabs by spaces in code fragments. # Minimum value: 1, maximum value: 16, default value: 4. -TAB_SIZE = 4 +TAB_SIZE = 2 # This tag can be used to specify a number of aliases that act as commands in # the documentation. An alias has the form: @@ -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 @@ -889,7 +890,7 @@ FILTER_SOURCE_PATTERNS = # (index.html). This can be useful if you have a project on for instance GitHub # and want to reuse the introduction page also for the doxygen output. -USE_MDFILE_AS_MAINPAGE = +USE_MDFILE_AS_MAINPAGE = doc/common/main_page.md #--------------------------------------------------------------------------- # Configuration options related to source browsing @@ -1147,7 +1148,7 @@ HTML_DYNAMIC_SECTIONS = NO # Minimum value: 0, maximum value: 9999, default value: 100. # This tag requires that the tag GENERATE_HTML is set to YES. -HTML_INDEX_NUM_ENTRIES = 100 +HTML_INDEX_NUM_ENTRIES = 2 # If the GENERATE_DOCSET tag is set to YES, additional index files will be # generated that can be used as input for Apple's Xcode 3 integrated development @@ -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 6bb42037..80d8dfb9 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -1344,7 +1344,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/Toplex_map/include/gudhi/Lazy_toplex_map.h b/src/Toplex_map/include/gudhi/Lazy_toplex_map.h index b0b3706e..d7bccdff 100644 --- a/src/Toplex_map/include/gudhi/Lazy_toplex_map.h +++ b/src/Toplex_map/include/gudhi/Lazy_toplex_map.h @@ -93,7 +93,7 @@ class Lazy_toplex_map { std::unordered_map<Vertex, std::size_t> gamma0_lbounds; std::unordered_map<Vertex, Simplex_ptr_set> t0; - bool empty_toplex; // Is the empty simplex a toplex ? + bool empty_toplex = true; // Is the empty simplex a toplex ? typedef boost::heap::fibonacci_heap<std::pair<std::size_t, Vertex>> PriorityQueue; PriorityQueue cleaning_priority; diff --git a/src/Witness_complex/example/example_nearest_landmark_table.cpp b/src/Witness_complex/example/example_nearest_landmark_table.cpp index acaf7c54..441900c1 100644 --- a/src/Witness_complex/example/example_nearest_landmark_table.cpp +++ b/src/Witness_complex/example/example_nearest_landmark_table.cpp @@ -19,22 +19,19 @@ int main(int argc, char * const argv[]) { using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>; Simplex_tree simplex_tree; - Nearest_landmark_table nlt; // Example contains 5 witnesses and 5 landmarks - Nearest_landmark_range w0 = {std::make_pair(0, 0), std::make_pair(1, 1), std::make_pair(2, 2), - std::make_pair(3, 3), std::make_pair(4, 4)}; nlt.push_back(w0); - Nearest_landmark_range w1 = {std::make_pair(1, 0), std::make_pair(2, 1), std::make_pair(3, 2), - std::make_pair(4, 3), std::make_pair(0, 4)}; nlt.push_back(w1); - Nearest_landmark_range w2 = {std::make_pair(2, 0), std::make_pair(3, 1), std::make_pair(4, 2), - std::make_pair(0, 3), std::make_pair(1, 4)}; nlt.push_back(w2); - Nearest_landmark_range w3 = {std::make_pair(3, 0), std::make_pair(4, 1), std::make_pair(0, 2), - std::make_pair(1, 3), std::make_pair(2, 4)}; nlt.push_back(w3); - Nearest_landmark_range w4 = {std::make_pair(4, 0), std::make_pair(0, 1), std::make_pair(1, 2), - std::make_pair(2, 3), std::make_pair(3, 4)}; nlt.push_back(w4); + Nearest_landmark_table nlt = { + {{0, 0.0}, {1, 0.1}, {2, 0.2}, {3, 0.3}, {4, 0.4}}, // witness 0 + {{1, 0.0}, {2, 0.1}, {3, 0.2}, {4, 0.3}, {0, 0.4}}, // witness 1 + {{2, 0.0}, {3, 0.1}, {4, 0.2}, {0, 0.3}, {1, 0.4}}, // witness 2 + {{3, 0.0}, {4, 0.1}, {0, 0.2}, {1, 0.3}, {2, 0.4}}, // witness 3 + {{4, 0.0}, {0, 0.1}, {1, 0.2}, {2, 0.3}, {3, 0.4}} // witness 4 + }; + /* distance(witness3, landmark3) is 0, distance(witness3, landmark4) is 0.1, etc. */ Witness_complex witness_complex(nlt); - witness_complex.create_complex(simplex_tree, 4.1); + witness_complex.create_complex(simplex_tree, .41); std::cout << "Number of simplices: " << simplex_tree.num_simplices() << std::endl; diff --git a/src/Witness_complex/include/gudhi/Strong_witness_complex.h b/src/Witness_complex/include/gudhi/Strong_witness_complex.h index fd6b3f38..03d6d2e7 100644 --- a/src/Witness_complex/include/gudhi/Strong_witness_complex.h +++ b/src/Witness_complex/include/gudhi/Strong_witness_complex.h @@ -40,10 +40,13 @@ namespace witness_complex { * \brief Constructs strong witness complex for a given table of nearest landmarks with respect to witnesses. * \ingroup witness_complex * - * \tparam Nearest_landmark_table_ needs to be a range of a range of pairs of nearest landmarks and distances. + * \tparam Nearest_landmark_table_ needs to be a range (one entry per witness) + * of sorted ranges of pairs of nearest landmarks and distances. * The class Nearest_landmark_table_::value_type must be a copiable range. * The range of pairs must admit a member type 'iterator'. The dereference type - * of the pair range iterator needs to be 'std::pair<std::size_t, double>'. + * of the pair range iterator needs to be 'std::pair<std::size_t, double>' + * where the first element is the index of the landmark, and the second its + * (squared) distance to the witness. */ template< class Nearest_landmark_table_ > class Strong_witness_complex { diff --git a/src/Witness_complex/include/gudhi/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h index 67885258..1f61f8f2 100644 --- a/src/Witness_complex/include/gudhi/Witness_complex.h +++ b/src/Witness_complex/include/gudhi/Witness_complex.h @@ -75,10 +75,13 @@ class Witness_complex { /** * \brief Initializes member variables before constructing simplicial complex. * \details Records nearest landmark table. - * @param[in] nearest_landmark_table needs to be a range of a range of pairs of nearest landmarks and distances. + * @param[in] nearest_landmark_table needs to be a range (one entry per witness) + * of sorted ranges of pairs of nearest landmarks and distances. * The class Nearest_landmark_table_::value_type must be a copiable range. * The range of pairs must admit a member type 'iterator'. The dereference type - * of the pair range iterator needs to be 'std::pair<std::size_t, double>'. + * of the pair range iterator needs to be 'std::pair<std::size_t, double>' + * where the first element is the index of the landmark, and the second its + * (squared) distance to the witness. */ Witness_complex(Nearest_landmark_table_ const & nearest_landmark_table) @@ -108,8 +111,8 @@ class Witness_complex { } ActiveWitnessList active_witnesses; Landmark_id k = 0; /* current dimension in iterative construction */ - for (auto w : nearest_landmark_table_) - active_witnesses.push_back(ActiveWitness(w)); + for (auto&& w : nearest_landmark_table_) + active_witnesses.emplace_back(w); while (!active_witnesses.empty() && k <= limit_dimension) { typename ActiveWitnessList::iterator aw_it = active_witnesses.begin(); std::vector<Landmark_id> simplex; 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/doc/main_page.h b/src/common/doc/main_page.md index afe9b68c..e61eee81 100644 --- a/src/common/doc/main_page.h +++ b/src/common/doc/main_page.md @@ -1,261 +1,289 @@ -/*! \mainpage The C++ library - * \tableofcontents - * \image html "Gudhi_banner.png" "" width=20cm - * - * \section Introduction Introduction - * The GUDHI library (Geometry Understanding in Higher Dimensions) is a generic open source - * <a class="el" target="_blank" href="http://gudhi.gforge.inria.fr/doc/latest/">C++ library</a> for - * Computational Topology and Topological Data Analysis - * (<a class="el" target="_blank" href="https://en.wikipedia.org/wiki/Topological_data_analysis">TDA</a>). - * The GUDHI library intends to help the development of new algorithmic solutions in TDA and their transfer to - * applications. It provides robust, efficient, flexible and easy to use implementations of state-of-the-art - * algorithms and data structures. - * - * The current release of the GUDHI library includes: - * - * \li Data structures to represent, construct and manipulate simplicial complexes. - * \li Simplification of simplicial complexes by edge contraction. - * \li Algorithms to compute persistent homology and bottleneck distance. - * - * All data-structures are generic and several of their aspects can be parameterized via template classes. - * We refer to \cite gudhilibrary_ICMS14 for a detailed description of the design of the library. - * - \section DataStructures Data structures - \subsection AlphaComplexDataStructure Alpha complex - \image html "alpha_complex_representation.png" "Alpha complex representation" -<table border="0"> +[TOC] + +# The C++ library {#main_page} +\image html "Gudhi_banner.png" +<br><br><br><br> + +## Complexes {#Complexes} +### Cubical complex + +<table> <tr> - <td width="25%"> - <b>Author:</b> Vincent Rouvreau<br> + <td width="35%" rowspan=2> + \image html "Cubical_complex_representation.png" + </td> + <td width="50%"> + The cubical complex is an example of a structured complex useful in computational mathematics (specially + rigorous numerics) and image analysis.<br> + </td> + <td width="15%"> + <b>Author:</b> Pawel Dlotko<br> <b>Introduced in:</b> GUDHI 1.3.0<br> <b>Copyright:</b> GPL v3<br> - <b>Requires:</b> \ref eigen3 and<br> - \ref cgal ≥ 4.7.0 for Alpha_complex<br> - \ref cgal ≥ 4.11.0 for Alpha_complex_3d </td> - <td width="75%"> - Alpha_complex is a simplicial complex constructed from the finite cells of a Delaunay Triangulation.<br> + </tr> + <tr> + <td colspan=2 height="25"> + <b>User manual:</b> \ref cubical_complex - <b>Reference manual:</b> Gudhi::cubical_complex::Bitmap_cubical_complex + </td> + </tr> +</table> + +### Simplicial complex + +#### Alpha complex + +<table> + <tr> + <td width="35%" rowspan=2> + \image html "alpha_complex_representation.png" + </td> + <td width="50%"> + Alpha complex is a simplicial complex constructed from the finite cells of a Delaunay Triangulation.<br> The filtration value of each simplex is computed as the square of the circumradius of the simplex if the circumsphere is empty (the simplex is then said to be Gabriel), and as the minimum of the filtration values of the codimension 1 cofaces that make it not Gabriel otherwise. All simplices that have a filtration value strictly greater than a given alpha squared value are not inserted into the complex.<br> + </td> + <td width="15%"> + <b>Author:</b> Vincent Rouvreau<br> + <b>Introduced in:</b> GUDHI 1.3.0<br> + <b>Copyright:</b> GPL v3<br> + <b>Requires:</b> \ref eigen3 and<br> + \ref cgal ≥ 4.7.0 for Alpha_complex<br> + \ref cgal ≥ 4.11.0 for Alpha_complex_3d + </td> + </tr> + <tr> + <td colspan=2 height="25"> <b>User manual:</b> \ref alpha_complex - <b>Reference manual:</b> Gudhi::alpha_complex::Alpha_complex and Gudhi::alpha_complex::Alpha_complex_3d </td> </tr> </table> - \subsection CechComplexDataStructure Čech complex - \image html "cech_complex_representation.png" "Čech complex representation" -<table border="0"> - <tr> - <td width="25%"> + +#### Čech complex + +<table> + <tr> + <td width="35%" rowspan=2> + \image html "cech_complex_representation.png" + </td> + <td width="50%"> + The Čech complex is a simplicial complex constructed from a proximity graph. + The set of all simplices is filtered by the radius of their minimal enclosing ball. + </td> + <td width="15%"> <b>Author:</b> Vincent Rouvreau<br> <b>Introduced in:</b> GUDHI 2.2.0<br> <b>Copyright:</b> GPL v3<br> </td> - <td width="75%"> - The Čech complex is a simplicial complex constructed from a proximity graph.<br> - The set of all simplices is filtered by the radius of their minimal enclosing ball.<br> + </tr> + <tr> + <td colspan=2 height="25"> <b>User manual:</b> \ref cech_complex - <b>Reference manual:</b> Gudhi::cech_complex::Cech_complex </td> </tr> </table> - \subsection CubicalComplexDataStructure Cubical complex - \image html "Cubical_complex_representation.png" "Cubical complex representation" -<table border="0"> + +#### Rips complex + +<table> <tr> - <td width="25%"> - <b>Author:</b> Pawel Dlotko<br> - <b>Introduced in:</b> GUDHI 1.3.0<br> - <b>Copyright:</b> GPL v3<br> + <td width="35%" rowspan=2> + \image html "rips_complex_representation.png" </td> - <td width="75%"> - The cubical complex is an example of a structured complex useful in computational mathematics (specially - rigorous numerics) and image analysis.<br> - <b>User manual:</b> \ref cubical_complex - <b>Reference manual:</b> Gudhi::cubical_complex::Bitmap_cubical_complex + <td width="50%"> + Rips complex is a simplicial complex constructed from a one skeleton graph.<br> + The filtration value of each edge is computed from a user-given distance function and is inserted until a + user-given threshold value.<br> + This complex can be built from a point cloud and a distance function, or from a distance matrix. </td> - </tr> -</table> - \subsection RipsComplexDataStructure Rips complex - \image html "rips_complex_representation.png" "Rips complex representation" -<table border="0"> - <tr> - <td width="25%"> + <td width="15%"> <b>Author:</b> Clément Maria, Pawel Dlotko, Vincent Rouvreau, Marc Glisse<br> <b>Introduced in:</b> GUDHI 2.0.0<br> <b>Copyright:</b> GPL v3<br> </td> - <td width="75%"> - Rips_complex is a simplicial complex constructed from a one skeleton graph.<br> - The filtration value of each edge is computed from a user-given distance function and is inserted until a - user-given threshold value.<br> - This complex can be built from a point cloud and a distance function, or from a distance matrix.<br> + </tr> + <tr> + <td colspan=2 height="25"> <b>User manual:</b> \ref rips_complex - <b>Reference manual:</b> Gudhi::rips_complex::Rips_complex </td> </tr> </table> - \subsection SimplexTreeDataStructure Simplex tree - \image html "Simplex_tree_representation.png" "Simplex tree representation" -<table border="0"> + +#### Witness complex + +<table> <tr> - <td width="25%"> - <b>Author:</b> Clément Maria<br> - <b>Introduced in:</b> GUDHI 1.0.0<br> + <td width="35%" rowspan=2> + \image html "Witness_complex_representation.png" + </td> + <td width="50%"> + Witness complex \f$ Wit(W,L) \f$ is a simplicial complex defined on two sets of points in \f$\mathbb{R}^D\f$. + The data structure is described in \cite boissonnatmariasimplextreealgorithmica . + </td> + <td width="15%"> + <b>Author:</b> Siargey Kachanovich<br> + <b>Introduced in:</b> GUDHI 1.3.0<br> <b>Copyright:</b> GPL v3<br> + <b>Euclidean version requires:</b> \ref cgal ≥ 4.6.0 and \ref eigen3 </td> - <td width="75%"> - The simplex tree is an efficient and flexible - data structure for representing general (filtered) simplicial complexes. The data structure - is described in \cite boissonnatmariasimplextreealgorithmica .<br> - <b>User manual:</b> \ref simplex_tree - <b>Reference manual:</b> Gudhi::Simplex_tree + </tr> + <tr> + <td colspan=2 height="25"> + <b>User manual:</b> \ref witness_complex - <b>Reference manual:</b> Gudhi::witness_complex::SimplicialComplexForWitness </td> </tr> </table> - \subsection CoverComplexDataStructure Cover Complexes - \image html "gicvisu.jpg" "Graph Induced Complex of a point cloud." -<table border="0"> + +### Cover Complexes +<table> <tr> - <td width="25%"> - <b>Author:</b> Mathieu Carrière<br> - <b>Introduced in:</b> GUDHI 2.1.0<br> - <b>Copyright:</b> GPL v3<br> - <b>Requires:</b> \ref cgal ≥ 4.8.1 + <td width="35%" rowspan=2> + \image html "gicvisu.jpg" </td> - <td width="75%"> + <td width="50%"> Nerves and Graph Induced Complexes are cover complexes, i.e. simplicial complexes that provably contain topological information about the input data. They can be computed with a cover of the data, that comes i.e. from the preimage of a family of intervals covering the image of a scalar-valued function defined on the data. <br> <b>User manual:</b> \ref cover_complex - <b>Reference manual:</b> Gudhi::cover_complex::Cover_complex </td> + <td width="15%"> + <b>Author:</b> Mathieu Carrière<br> + <b>Introduced in:</b> GUDHI 2.1.0<br> + <b>Copyright:</b> GPL v3<br> + <b>Requires:</b> \ref cgal ≥ 4.8.1 + </td> + </tr> + <tr> + <td colspan=2 height="25"> + <b>User manual:</b> \ref cover_complex - <b>Reference manual:</b> Gudhi::cover_complex::Cover_complex + </td> </tr> </table> - \subsection SkeletonBlockerDataStructure Skeleton blocker - \image html "ds_representation.png" "Skeleton blocker representation" -<table border="0"> + +## Data structures and basic operations {#DataStructuresAndBasicOperations} + +### Data structures + +#### Simplex tree +<table> <tr> - <td width="25%"> - <b>Author:</b> David Salinas<br> - <b>Introduced in:</b> GUDHI 1.1.0<br> + <td width="35%" rowspan=2> + \image html "Simplex_tree_representation.png" + </td> + <td width="50%"> + The simplex tree is an efficient and flexible + data structure for representing general (filtered) simplicial complexes. The data structure + is described in \cite boissonnatmariasimplextreealgorithmica . + </td> + <td width="15%"> + <b>Author:</b> Clément Maria<br> + <b>Introduced in:</b> GUDHI 1.0.0<br> <b>Copyright:</b> GPL v3<br> </td> - <td width="75%"> + </tr> + <tr> + <td colspan=2 height="25"> + <b>User manual:</b> \ref simplex_tree - <b>Reference manual:</b> Gudhi::Simplex_tree + </td> + </tr> +</table> + +#### Skeleton blocker + +<table> + <tr> + <td width="35%" rowspan=2> + \image html "ds_representation.png" + </td> + <td width="50%"> The Skeleton-Blocker data-structure proposes a light encoding for simplicial complexes by storing only an *implicit* representation of its simplices \cite socg_blockers_2011,\cite blockers2012. Intuitively, it just stores the 1-skeleton of a simplicial complex with a graph and the set of its "missing faces" that is very small in practice. This data-structure handles all simplicial complexes operations such as simplex enumeration or simplex removal but operations that are particularly efficient are operations that do not require simplex enumeration such as edge - iteration, link computation or simplex contraction.<br> - <b>User manual:</b> \ref skbl - <b>Reference manual:</b> Gudhi::skeleton_blocker::Skeleton_blocker_complex + iteration, link computation or simplex contraction. </td> - </tr> -</table> - \subsection TangentialComplexDataStructure Tangential complex - \image html "tc_examples.png" "Tangential complex representation" -<table border="0"> - <tr> - <td width="25%"> - <b>Author:</b> Clément Jamin<br> - <b>Introduced in:</b> GUDHI 2.0.0<br> + <td width="15%"> + <b>Author:</b> David Salinas<br> + <b>Introduced in:</b> GUDHI 1.1.0<br> <b>Copyright:</b> GPL v3<br> - <b>Requires:</b> \ref cgal ≥ 4.8.1 and \ref eigen3 </td> - <td width="75%"> - A Tangential Delaunay complex is a <a target="_blank" href="https://en.wikipedia.org/wiki/Simplicial_complex">simplicial complex</a> - designed to reconstruct a \f$ k \f$-dimensional manifold embedded in \f$ d \f$-dimensional Euclidean space. - The input is a point sample coming from an unknown manifold. - The running time depends only linearly on the extrinsic dimension \f$ d \f$ - and exponentially on the intrinsic dimension \f$ k \f$.<br> - <b>User manual:</b> \ref tangential_complex - <b>Reference manual:</b> Gudhi::tangential_complex::Tangential_complex + </tr> + <tr> + <td colspan=2 height="25"> + <b>User manual:</b> \ref skbl - <b>Reference manual:</b> Gudhi::skeleton_blocker::Skeleton_blocker_complex </td> </tr> </table> - \subsection ToplexMapDataStructure Toplex Map - \image html "map.png" "Toplex map representation" -<table border="0"> + +#### Toplex Map + +<table> <tr> - <td width="25%"> - <b>Author:</b> François Godi<br> - <b>Introduced in:</b> GUDHI 2.1.0<br> - <b>Copyright:</b> GPL v3<br> + <td width="35%" rowspan=2> + \image html "map.png" </td> - <td width="75%"> + <td width="50%"> The Toplex map data structure is composed firstly of a raw storage of toplices (the maximal simplices) and secondly of a map which associate any vertex to a set of pointers toward all toplices containing this vertex. - <b>User manual:</b> \ref toplex_map - <b>Reference manual:</b> Gudhi::Toplex_map </td> - </tr> - - \subsection WitnessComplexDataStructure Witness complex - \image html "Witness_complex_representation.png" "Witness complex representation" -<table border="0"> - <tr> - <td width="25%"> - <b>Author:</b> Siargey Kachanovich<br> - <b>Introduced in:</b> GUDHI 1.3.0<br> - <b>Copyright:</b> GPL v3<br> - <b>Euclidean version requires:</b> \ref cgal ≥ 4.6.0 and \ref eigen3 - </td> - <td width="75%"> - Witness complex \f$ Wit(W,L) \f$ is a simplicial complex defined on two sets of points in \f$\mathbb{R}^D\f$. - The data structure is described in \cite boissonnatmariasimplextreealgorithmica .<br> - <b>User manual:</b> \ref witness_complex - <b>Reference manual:</b> Gudhi::witness_complex::SimplicialComplexForWitness - </td> - </tr> -</table> - - \section Toolbox Toolbox - - \subsection BottleneckDistanceToolbox Bottleneck distance - \image html "perturb_pd.png" "Bottleneck distance is the length of the longest edge" -<table border="0"> - <tr> - <td width="25%"> + <td width="15%"> <b>Author:</b> François Godi<br> - <b>Introduced in:</b> GUDHI 2.0.0<br> + <b>Introduced in:</b> GUDHI 2.1.0<br> <b>Copyright:</b> GPL v3<br> - <b>Requires:</b> \ref cgal ≥ 4.8.1 </td> - <td width="75%"> - Bottleneck distance measures the similarity between two persistence diagrams. - It's the shortest distance b for which there exists a perfect matching between - the points of the two diagrams (+ all the diagonal points) such that - any couple of matched points are at distance at most b. - <br> - <b>User manual:</b> \ref bottleneck_distance + </tr> + <tr> + <td colspan=2 height="25"> + <b>User manual:</b> \ref toplex_map - <b>Reference manual:</b> Gudhi::Toplex_map </td> </tr> </table> - \subsection ContractionToolbox Contraction - \image html "sphere_contraction_representation.png" "Sphere contraction example" -<table border="0"> + +### Basic operations + +#### Contraction + +<table> <tr> - <td width="25%"> + <td width="35%" rowspan=2> + \image html "sphere_contraction_representation.png" + </td> + <td width="15%"> <b>Author:</b> David Salinas<br> <b>Introduced in:</b> GUDHI 1.1.0<br> <b>Copyright:</b> GPL v3<br> </td> - <td width="75%"> + <td width="50%"> The purpose of this package is to offer a user-friendly interface for edge contraction simplification of huge simplicial complexes. It uses the \ref skbl data-structure whose size remains small during simplification of most used geometrical complexes of topological data analysis such as the Rips or the Delaunay complexes. In practice, - the size of this data-structure is even much lower than the total number of simplices.<br> + the size of this data-structure is even much lower than the total number of simplices. + </td> + </tr> + <tr> + <td colspan=2 height="25"> <b>User manual:</b> \ref contr </td> </tr> </table> - \subsection PersistentCohomologyToolbox Persistent Cohomology - \image html "3DTorus_poch.png" "Rips Persistent Cohomology on a 3D Torus" -<table border="0"> + +## Topological descriptors computation {#TopologicalDescriptorsComputation} + +### Persistent Cohomology + +<table> <tr> - <td width="25%"> - <b>Author:</b> Clément Maria<br> - <b>Introduced in:</b> GUDHI 1.0.0<br> - <b>Copyright:</b> GPL v3<br> + <td width="35%" rowspan=2> + \image html "3DTorus_poch.png" </td> - <td width="75%"> + <td width="50%"> The theory of homology consists in attaching to a topological space a sequence of (homology) groups, capturing global topological features like connected components, holes, cavities, etc. Persistent homology studies the evolution -- birth, life and death -- of these features when the topological space is changing. Consequently, the @@ -263,27 +291,101 @@ scheme. Computation of persistent cohomology using the algorithm of \cite DBLP:journals/dcg/SilvaMV11 and \cite DBLP:journals/corr/abs-1208-5018 and the Compressed Annotation Matrix implementation of - \cite DBLP:conf/esa/BoissonnatDM13 .<br> + \cite DBLP:conf/esa/BoissonnatDM13 . + </td> + <td width="15%"> + <b>Author:</b> Clément Maria<br> + <b>Introduced in:</b> GUDHI 1.0.0<br> + <b>Copyright:</b> GPL v3<br> + </td> + </tr> + <tr> + <td colspan=2 height="25"> <b>User manual:</b> \ref persistent_cohomology - <b>Reference manual:</b> Gudhi::persistent_cohomology::Persistent_cohomology </td> </tr> -</table> - \subsection PersistenceRepresentationsToolbox Persistence representations - \image html "average_landscape.png" "Persistence representations" -<table border="0"> +</table> + +## Manifold reconstruction {#ManifoldReconstruction} + +### Tangential complex + +<table> + <tr> + <td width="35%" rowspan=2> + \image html "tc_examples.png" + </td> + <td width="50%"> + A Tangential Delaunay complex is a <a target="_blank" href="https://en.wikipedia.org/wiki/Simplicial_complex">simplicial complex</a> + designed to reconstruct a \f$ k \f$-dimensional manifold embedded in \f$ d \f$-dimensional Euclidean space. + The input is a point sample coming from an unknown manifold. + The running time depends only linearly on the extrinsic dimension \f$ d \f$ + and exponentially on the intrinsic dimension \f$ k \f$. + </td> + <td width="15%"> + <b>Author:</b> Clément Jamin<br> + <b>Introduced in:</b> GUDHI 2.0.0<br> + <b>Copyright:</b> GPL v3<br> + <b>Requires:</b> \ref cgal ≥ 4.8.1 and \ref eigen3 + </td> + </tr> + <tr> + <td colspan=2 height="25"> + <b>User manual:</b> \ref tangential_complex - <b>Reference manual:</b> Gudhi::tangential_complex::Tangential_complex + </td> + </tr> +</table> + +## Topological descriptors tools {#TopologicalDescriptorsTools} + +### Bottleneck distance + +<table> + <tr> + <td width="35%" rowspan=2> + \image html "perturb_pd.png" + </td> + <td width="50%"> + Bottleneck distance measures the similarity between two persistence diagrams. + It's the shortest distance b for which there exists a perfect matching between + the points of the two diagrams (+ all the diagonal points) such that + any couple of matched points are at distance at most b. + </td> + <td width="15%"> + <b>Author:</b> François Godi<br> + <b>Introduced in:</b> GUDHI 2.0.0<br> + <b>Copyright:</b> GPL v3<br> + <b>Requires:</b> \ref cgal ≥ 4.8.1 + </td> + </tr> + <tr> + <td colspan=2 height="25"> + <b>User manual:</b> \ref bottleneck_distance + </td> + </tr> +</table> + +### Persistence representations + +<table> <tr> - <td width="25%"> + <td width="35%" rowspan=2> + \image html "average_landscape.png" + </td> + <td width="50%"> + It contains implementation of various representations of persistence diagrams; diagrams themselves, persistence + landscapes (rigorous and grid version), persistence heath maps, vectors and others. It implements basic + functionalities which are neccessary to use persistence in statistics and machine learning. + </td> + <td width="15%"> <b>Author:</b> Pawel Dlotko<br> <b>Introduced in:</b> GUDHI 2.1.0<br> <b>Copyright:</b> GPL v3<br> </td> - <td width="75%"> - It contains implementation of various representations of persistence diagrams; diagrams themselves, persistence - landscapes (rigorous and grid version), persistence heath maps, vectors and others. It implements basic - functionalities which are neccessary to use persistence in statistics and machine learning.<br> + </tr> + <tr> + <td colspan=2 height="25"> <b>User manual:</b> \ref Persistence_representations </td> </tr> </table> - -*/ 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/cython/strong_witness_complex.pyx b/src/cython/cython/strong_witness_complex.pyx index 74c5cb05..4b7bff34 100644 --- a/src/cython/cython/strong_witness_complex.pyx +++ b/src/cython/cython/strong_witness_complex.pyx @@ -47,8 +47,10 @@ cdef class StrongWitnessComplex: def __init__(self, nearest_landmark_table=None): """StrongWitnessComplex constructor. - :param nearest_landmark_table: A list of nearest landmark. - :type nearest_landmark_table: list of list of pair of unsigned and double + :param nearest_landmark_table: A list of lists of nearest landmarks and their distances. + `nearest_landmark_table[w][k]==(l,d)` means that l is the k-th nearest landmark to + witness w, and d is the (squared) distance between l and w. + :type nearest_landmark_table: list of list of pair of int and float """ # The real cython constructor @@ -65,10 +67,10 @@ cdef class StrongWitnessComplex: """ return self.thisptr != NULL - def create_simplex_tree(self, max_alpha_square, limit_dimension = -1): + def create_simplex_tree(self, max_alpha_square = float('inf'), limit_dimension = -1): """ - :param max_alpha_square: The maximum alpha square threshold the - simplices shall not exceed. Default is set to infinity. + :param max_alpha_square: The maximum relaxation parameter. + Default is set to infinity. :type max_alpha_square: float :returns: A simplex tree created from the Delaunay Triangulation. :rtype: SimplexTree diff --git a/src/cython/cython/witness_complex.pyx b/src/cython/cython/witness_complex.pyx index 8591465a..b1cce83f 100644 --- a/src/cython/cython/witness_complex.pyx +++ b/src/cython/cython/witness_complex.pyx @@ -47,8 +47,10 @@ cdef class WitnessComplex: def __init__(self, nearest_landmark_table=None): """WitnessComplex constructor. - :param nearest_landmark_table: A list of nearest landmark. - :type nearest_landmark_table: list of list of pair of unsigned and double + :param nearest_landmark_table: A list of lists of nearest landmarks and their distances. + `nearest_landmark_table[w][k]==(l,d)` means that l is the k-th nearest landmark to + witness w, and d is the (squared) distance between l and w. + :type nearest_landmark_table: list of list of pair of int and float """ # The real cython constructor @@ -65,10 +67,10 @@ cdef class WitnessComplex: """ return self.thisptr != NULL - def create_simplex_tree(self, max_alpha_square, limit_dimension = -1): + def create_simplex_tree(self, max_alpha_square = float('inf'), limit_dimension = -1): """ - :param max_alpha_square: The maximum alpha square threshold the - simplices shall not exceed. Default is set to infinity. + :param max_alpha_square: The maximum relaxation parameter. + Default is set to infinity. :type max_alpha_square: float :returns: A simplex tree created from the Delaunay Triangulation. :rtype: SimplexTree diff --git a/src/cython/doc/alpha_complex_sum.inc b/src/cython/doc/alpha_complex_sum.inc index 1680a712..806988bb 100644 --- a/src/cython/doc/alpha_complex_sum.inc +++ b/src/cython/doc/alpha_complex_sum.inc @@ -1,22 +1,20 @@ -================================================================= =================================== =================================== -:Author: Vincent Rouvreau :Introduced in: GUDHI 2.0.0 :Copyright: GPL v3 -:Requires: CGAL :math:`\geq` 4.7.0 Eigen3 -================================================================= =================================== =================================== +.. table:: + :widths: 30 50 20 -+----------------------------------------------------------------+------------------------------------------------------------------------+ -| .. figure:: | Alpha_complex is a simplicial complex constructed from the finite | -| ../../doc/Alpha_complex/alpha_complex_representation.png | cells of a Delaunay Triangulation. | -| :alt: Alpha complex representation | | -| :figclass: align-center | The filtration value of each simplex is computed as the square of the | -| | circumradius of the simplex if the circumsphere is empty (the simplex | -| Alpha complex representation | is then said to be Gabriel), and as the minimum of the filtration | -| | values of the codimension 1 cofaces that make it not Gabriel | -| | otherwise. All simplices that have a filtration value strictly | -| | greater than a given alpha squared value are not inserted into the | -| | complex. | -| | | -| | This package requires having CGAL version 4.7 or higher (4.8.1 is | -| | advised for better performance). | -+----------------------------------------------------------------+------------------------------------------------------------------------+ -| :doc:`alpha_complex_user` | :doc:`alpha_complex_ref` | -+----------------------------------------------------------------+------------------------------------------------------------------------+ + +----------------------------------------------------------------+------------------------------------------------------------------------+-----------------------------------------------+ + | .. figure:: | Alpha complex is a simplicial complex constructed from the finite | :Author: Vincent Rouvreau | + | ../../doc/Alpha_complex/alpha_complex_representation.png | cells of a Delaunay Triangulation. | | + | :alt: Alpha complex representation | | :Introduced in: GUDHI 2.0.0 | + | :figclass: align-center | The filtration value of each simplex is computed as the square of the | | + | | circumradius of the simplex if the circumsphere is empty (the simplex | :Copyright: GPL v3 | + | | is then said to be Gabriel), and as the minimum of the filtration | | + | | values of the codimension 1 cofaces that make it not Gabriel | :Requires: Eigen3 and CGAL :math:`\geq` 4.7.0 | + | | otherwise. All simplices that have a filtration value strictly | | + | | greater than a given alpha squared value are not inserted into the | | + | | complex. | | + | | | | + | | This package requires having CGAL version 4.7 or higher (4.8.1 is | | + | | advised for better performance). | | + +----------------------------------------------------------------+------------------------------------------------------------------------+-----------------------------------------------+ + | * :doc:`alpha_complex_user` | * :doc:`alpha_complex_ref` | + +----------------------------------------------------------------+------------------------------------------------------------------------------------------------------------------------+ diff --git a/src/cython/doc/bottleneck_distance_sum.inc b/src/cython/doc/bottleneck_distance_sum.inc index 030fad9e..41b9c5a3 100644 --- a/src/cython/doc/bottleneck_distance_sum.inc +++ b/src/cython/doc/bottleneck_distance_sum.inc @@ -1,15 +1,14 @@ -================================================================= =================================== =================================== -:Author: François Godi :Introduced in: GUDHI 2.0.0 :Copyright: GPL v3 -:Requires: CGAL :math:`\geq` 4.8.0 -================================================================= =================================== =================================== +.. table:: + :widths: 30 50 20 -+-----------------------------------------------------------------+----------------------------------------------------------------------+ -| .. figure:: | Bottleneck distance measures the similarity between two persistence | -| ../../doc/Bottleneck_distance/perturb_pd.png | diagrams. It's the shortest distance b for which there exists a | -| :figclass: align-center | perfect matching between the points of the two diagrams (+ all the | -| | diagonal points) such that any couple of matched points are at | -| Bottleneck distance is the length of | distance at most b. | -| the longest edge | | -+-----------------------------------------------------------------+----------------------------------------------------------------------+ -| :doc:`bottleneck_distance_user` | | -+-----------------------------------------------------------------+----------------------------------------------------------------------+ + +-----------------------------------------------------------------+----------------------------------------------------------------------+-----------------------------------------------+ + | .. figure:: | Bottleneck distance measures the similarity between two persistence | :Author: François Godi | + | ../../doc/Bottleneck_distance/perturb_pd.png | diagrams. It's the shortest distance b for which there exists a | | + | :figclass: align-center | perfect matching between the points of the two diagrams (+ all the | :Introduced in: GUDHI 2.0.0 | + | | diagonal points) such that any couple of matched points are at | | + | Bottleneck distance is the length of | distance at most b. | :Copyright: GPL v3 | + | the longest edge | | | + | | | :Requires: CGAL :math:`\geq` 4.8.0 | + +-----------------------------------------------------------------+----------------------------------------------------------------------+-----------------------------------------------+ + | * :doc:`bottleneck_distance_user` | | + +-----------------------------------------------------------------+----------------------------------------------------------------------------------------------------------------------+ diff --git a/src/cython/doc/conf.py b/src/cython/doc/conf.py index 4a54d4fd..ce08f679 100755 --- a/src/cython/doc/conf.py +++ b/src/cython/doc/conf.py @@ -125,7 +125,7 @@ html_theme_options = { "sidebarbgcolor": "#A1ADCD", "sidebartextcolor": "black", "sidebarlinkcolor": "#334D5C", - "body_max_width": "1200px", + "body_max_width": "100%", } # Add any paths that contain custom themes here, relative to this directory. diff --git a/src/cython/doc/cubical_complex_sum.inc b/src/cython/doc/cubical_complex_sum.inc index 280ad0e0..6dcf8e48 100644 --- a/src/cython/doc/cubical_complex_sum.inc +++ b/src/cython/doc/cubical_complex_sum.inc @@ -1,15 +1,14 @@ -================================================================= =================================== =================================== -:Author: Pawel Dlotko :Introduced in: GUDHI 2.0.0 :Copyright: GPL v3 -================================================================= =================================== =================================== +.. table:: + :widths: 30 50 20 -+--------------------------------------------------------------------------+----------------------------------------------------------------------+ -| .. figure:: | The cubical complex is an example of a structured complex useful in | -| ../../doc/Bitmap_cubical_complex/Cubical_complex_representation.png | computational mathematics (specially rigorous numerics) and image | -| :alt: Cubical complex representation | analysis. | -| :figclass: align-center | | -| | | -| Cubical complex representation | | -+--------------------------------------------------------------------------+----------------------------------------------------------------------+ -| :doc:`cubical_complex_user` | * :doc:`cubical_complex_ref` | -| | * :doc:`periodic_cubical_complex_ref` | -+--------------------------------------------------------------------------+----------------------------------------------------------------------+ + +--------------------------------------------------------------------------+----------------------------------------------------------------------+-----------------------------+ + | .. figure:: | The cubical complex is an example of a structured complex useful in | :Author: Pawel Dlotko | + | ../../doc/Bitmap_cubical_complex/Cubical_complex_representation.png | computational mathematics (specially rigorous numerics) and image | | + | :alt: Cubical complex representation | analysis. | :Introduced in: GUDHI 2.0.0 | + | :figclass: align-center | | | + | | | :Copyright: GPL v3 | + | | | | + +--------------------------------------------------------------------------+----------------------------------------------------------------------+-----------------------------+ + | * :doc:`cubical_complex_user` | * :doc:`cubical_complex_ref` | + | | * :doc:`periodic_cubical_complex_ref` | + +--------------------------------------------------------------------------+----------------------------------------------------------------------------------------------------+ 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/index.rst b/src/cython/doc/index.rst index 15cbe267..e379bc23 100644 --- a/src/cython/doc/index.rst +++ b/src/cython/doc/index.rst @@ -6,80 +6,73 @@ GUDHI Python module documentation :alt: Gudhi banner :figclass: align-center -Introduction -************ - -The Python interface for the Gudhi library (Geometry Understanding in Higher -Dimensions) is a generic open source -`Python module <http://gudhi.gforge.inria.fr/python/latest/>`_, for -Computational Topology and Topological Data Analysis -(`TDA <https://en.wikipedia.org/wiki/Topological_data_analysis>`_). -The GUDHI library intends to help the development of new algorithmic solutions -in TDA and their transfer to applications. It provides robust, efficient, -flexible and easy to use implementations of state-of-the-art algorithms and -data structures. - -The current release of the GUDHI library includes: +Complexes +********* -* Data structures to represent, construct and manipulate simplicial complexes. -* Simplification of simplicial complexes by edge contraction. -* Algorithms to compute persistent homology and bottleneck distance. +Cubical complexes +================= -We refer to :cite:`gudhilibrary_ICMS14` for a detailed description of the -design of the library. +.. include:: cubical_complex_sum.inc -Data structures -*************** +Simplicial complexes +==================== Alpha complex -============= +------------- .. include:: alpha_complex_sum.inc -Cover complexes -=============== +Rips complex +------------- + +.. include:: rips_complex_sum.inc -.. include:: nerve_gic_complex_sum.rst +Witness complex +--------------- -Cubical complex +.. include:: witness_complex_sum.inc + +Cover complexes =============== -.. include:: cubical_complex_sum.inc +.. include:: nerve_gic_complex_sum.inc -Rips complex -============ +Data structures and basic operations +************************************ -.. include:: rips_complex_sum.inc +Data structures +=============== Simplex tree -============ +------------ .. include:: simplex_tree_sum.inc +Topological descriptors computation +*********************************** + +Persistence cohomology +====================== + +.. include:: persistent_cohomology_sum.inc + +Manifold reconstruction +*********************** + Tangential complex ================== .. include:: tangential_complex_sum.inc -Witness complex -=============== - -.. include:: witness_complex_sum.inc - -Toolbox -******* +Topological descriptors tools +***************************** Bottleneck distance =================== .. include:: bottleneck_distance_sum.inc -Persistence cohomology -====================== - -.. include:: persistent_cohomology_sum.inc - Persistence graphical tools =========================== diff --git a/src/cython/doc/nerve_gic_complex_sum.inc b/src/cython/doc/nerve_gic_complex_sum.inc new file mode 100644 index 00000000..0e606fe1 --- /dev/null +++ b/src/cython/doc/nerve_gic_complex_sum.inc @@ -0,0 +1,16 @@ +.. table:: + :widths: 30 50 20 + + +----------------------------------------------------------------+------------------------------------------------------------------------+------------------------------------+ + | .. figure:: | Nerves and Graph Induced Complexes are cover complexes, i.e. | :Author: Mathieu Carrière | + | ../../doc/Nerve_GIC/gicvisu.jpg | simplicial complexes that provably contain topological information | | + | :alt: Graph Induced Complex of a point cloud. | about the input data. They can be computed with a cover of the data, | :Introduced in: GUDHI 2.3.0 | + | :figclass: align-center | that comes i.e. from the preimage of a family of intervals covering | | + | | the image of a scalar-valued function defined on the data. | :Copyright: GPL v3 | + | | | | + | | | :Requires: CGAL :math:`\geq` 4.8.1 | + | | | | + | | | | + +----------------------------------------------------------------+------------------------------------------------------------------------+------------------------------------+ + | * :doc:`nerve_gic_complex_user` | * :doc:`nerve_gic_complex_ref` | + +----------------------------------------------------------------+-------------------------------------------------------------------------------------------------------------+ diff --git a/src/cython/doc/nerve_gic_complex_sum.rst b/src/cython/doc/nerve_gic_complex_sum.rst deleted file mode 100644 index 523c119f..00000000 --- a/src/cython/doc/nerve_gic_complex_sum.rst +++ /dev/null @@ -1,15 +0,0 @@ -================================================================= =================================== =================================== -:Author: Mathieu Carrière :Introduced in: GUDHI 2.3.0 :Copyright: GPL v3 -:Requires: CGAL :math:`\geq` 4.8.1 -================================================================= =================================== =================================== - -+----------------------------------------------------------------+------------------------------------------------------------------------+ -| .. figure:: | Nerves and Graph Induced Complexes are cover complexes, i.e. | -| ../../doc/Nerve_GIC/gicvisu.jpg | simplicial complexes that provably contain topological information | -| :alt: Graph Induced Complex of a point cloud. | about the input data. They can be computed with a cover of the data, | -| :figclass: align-center | that comes i.e. from the preimage of a family of intervals covering | -| | the image of a scalar-valued function defined on the data. | -| Graph Induced Complex of a point cloud. | | -+----------------------------------------------------------------+------------------------------------------------------------------------+ -| :doc:`nerve_gic_complex_user` | :doc:`nerve_gic_complex_ref` | -+----------------------------------------------------------------+------------------------------------------------------------------------+ diff --git a/src/cython/doc/nerve_gic_complex_user.rst b/src/cython/doc/nerve_gic_complex_user.rst index 44f30e1a..9101f45d 100644 --- a/src/cython/doc/nerve_gic_complex_user.rst +++ b/src/cython/doc/nerve_gic_complex_user.rst @@ -7,14 +7,13 @@ Cover complexes user manual Definition ---------- -.. include:: nerve_gic_complex_sum.rst +.. include:: nerve_gic_complex_sum.inc 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/doc/persistence_graphical_tools_sum.inc b/src/cython/doc/persistence_graphical_tools_sum.inc index 5577cf99..b412de56 100644 --- a/src/cython/doc/persistence_graphical_tools_sum.inc +++ b/src/cython/doc/persistence_graphical_tools_sum.inc @@ -1,12 +1,14 @@ -================================================================= =================================== =================================== -:Author: Vincent Rouvreau :Introduced in: GUDHI 2.0.0 :Copyright: GPL v3 -:Requires: matplotlib numpy scipy -================================================================= =================================== =================================== +.. table:: + :widths: 30 50 20 -+-----------------------------------------------------------------+-----------------------------------------------------------------------+ -| .. figure:: | These graphical tools comes on top of persistence results and allows | -| img/graphical_tools_representation.png | the user to build easily persistence barcode, diagram or density. | -| | | -+-----------------------------------------------------------------+-----------------------------------------------------------------------+ -| :doc:`persistence_graphical_tools_user` | :doc:`persistence_graphical_tools_ref` | -+-----------------------------------------------------------------+-----------------------------------------------------------------------+ + +-----------------------------------------------------------------+-----------------------------------------------------------------------+-----------------------------------------------+ + | .. figure:: | These graphical tools comes on top of persistence results and allows | :Author: Vincent Rouvreau | + | img/graphical_tools_representation.png | the user to build easily persistence barcode, diagram or density. | | + | | | :Introduced in: GUDHI 2.0.0 | + | | | | + | | | :Copyright: GPL v3 | + | | | | + | | | :Requires: matplotlib, numpy and scipy | + +-----------------------------------------------------------------+-----------------------------------------------------------------------+-----------------------------------------------+ + | * :doc:`persistence_graphical_tools_user` | * :doc:`persistence_graphical_tools_ref` | + +-----------------------------------------------------------------+-----------------------------------------------------------------------------------------------------------------------+ diff --git a/src/cython/doc/persistent_cohomology_sum.inc b/src/cython/doc/persistent_cohomology_sum.inc index a26df1dc..20ca073c 100644 --- a/src/cython/doc/persistent_cohomology_sum.inc +++ b/src/cython/doc/persistent_cohomology_sum.inc @@ -1,27 +1,26 @@ -================================================================= =================================== =================================== -:Author: Clément Maria :Introduced in: GUDHI 2.0.0 :Copyright: GPL v3 -================================================================= =================================== =================================== +.. table:: + :widths: 30 50 20 -+-----------------------------------------------------------------+-----------------------------------------------------------------------+ -| .. figure:: | The theory of homology consists in attaching to a topological space | -| ../../doc/Persistent_cohomology/3DTorus_poch.png | a sequence of (homology) groups, capturing global topological | -| :figclass: align-center | features like connected components, holes, cavities, etc. Persistent | -| | homology studies the evolution -- birth, life and death -- of these | -| Rips Persistent Cohomology on a 3D | features when the topological space is changing. Consequently, the | -| Torus | theory is essentially composed of three elements: topological spaces, | -| | their homology groups and an evolution scheme. | -| | | -| | Computation of persistent cohomology using the algorithm of | -| | :cite:`DBLP:journals/dcg/SilvaMV11` and | -| | :cite:`DBLP:journals/corr/abs-1208-5018` and the Compressed | -| | Annotation Matrix implementation of | -| | :cite:`DBLP:conf/esa/BoissonnatDM13`. | -| | | -+-----------------------------------------------------------------+-----------------------------------------------------------------------+ -| :doc:`persistent_cohomology_user` | Please refer to each data structure that contains persistence | -| | feature for reference: | -| | | -| | * :doc:`simplex_tree_ref` | -| | * :doc:`cubical_complex_ref` | -| | * :doc:`periodic_cubical_complex_ref` | -+-----------------------------------------------------------------+-----------------------------------------------------------------------+ + +-----------------------------------------------------------------+-----------------------------------------------------------------------+-----------------------------------------------+ + | .. figure:: | The theory of homology consists in attaching to a topological space | :Author: Clément Maria | + | ../../doc/Persistent_cohomology/3DTorus_poch.png | a sequence of (homology) groups, capturing global topological | | + | :figclass: align-center | features like connected components, holes, cavities, etc. Persistent | :Introduced in: GUDHI 2.0.0 | + | | homology studies the evolution -- birth, life and death -- of these | | + | Rips Persistent Cohomology on a 3D | features when the topological space is changing. Consequently, the | :Copyright: GPL v3 | + | Torus | theory is essentially composed of three elements: topological spaces, | | + | | their homology groups and an evolution scheme. | | + | | | | + | | Computation of persistent cohomology using the algorithm of | | + | | :cite:`DBLP:journals/dcg/SilvaMV11` and | | + | | :cite:`DBLP:journals/corr/abs-1208-5018` and the Compressed | | + | | Annotation Matrix implementation of | | + | | :cite:`DBLP:conf/esa/BoissonnatDM13`. | | + | | | | + +-----------------------------------------------------------------+-----------------------------------------------------------------------+-----------------------------------------------+ + | * :doc:`persistent_cohomology_user` | Please refer to each data structure that contains persistence | + | | feature for reference: | + | | | + | | * :doc:`simplex_tree_ref` | + | | * :doc:`cubical_complex_ref` | + | | * :doc:`periodic_cubical_complex_ref` | + +-----------------------------------------------------------------+-----------------------------------------------------------------------------------------------------------------------+ diff --git a/src/cython/doc/rips_complex_sum.inc b/src/cython/doc/rips_complex_sum.inc index ea26769a..e8e505e2 100644 --- a/src/cython/doc/rips_complex_sum.inc +++ b/src/cython/doc/rips_complex_sum.inc @@ -1,17 +1,16 @@ -===================================================================== =========================== =================================== -:Author: Clément Maria, Pawel Dlotko, Vincent Rouvreau, Marc Glisse :Introduced in: GUDHI 2.0.0 :Copyright: GPL v3 -===================================================================== =========================== =================================== +.. table:: + :widths: 30 50 20 -+----------------------------------------------------------------+------------------------------------------------------------------------+ -| .. figure:: | Rips complex is a simplicial complex constructed from a one skeleton | -| ../../doc/Rips_complex/rips_complex_representation.png | graph. | -| :figclass: align-center | | -| | The filtration value of each edge is computed from a user-given | -| Rips complex representation | distance function and is inserted until a user-given threshold | -| | value. | -| | | -| | This complex can be built from a point cloud and a distance function, | -| | or from a distance matrix. | -+----------------------------------------------------------------+------------------------------------------------------------------------+ -| :doc:`rips_complex_user` | :doc:`rips_complex_ref` | -+----------------------------------------------------------------+------------------------------------------------------------------------+ + +----------------------------------------------------------------+------------------------------------------------------------------------+----------------------------------------------------------------------+ + | .. figure:: | Rips complex is a simplicial complex constructed from a one skeleton | :Authors: Clément Maria, Pawel Dlotko, Vincent Rouvreau, Marc Glisse | + | ../../doc/Rips_complex/rips_complex_representation.png | graph. | | + | :figclass: align-center | | :Introduced in: GUDHI 2.0.0 | + | | The filtration value of each edge is computed from a user-given | | + | | distance function and is inserted until a user-given threshold | :Copyright: GPL v3 | + | | value. | | + | | | | + | | This complex can be built from a point cloud and a distance function, | | + | | or from a distance matrix. | | + +----------------------------------------------------------------+------------------------------------------------------------------------+----------------------------------------------------------------------+ + | * :doc:`rips_complex_user` | * :doc:`rips_complex_ref` | + +----------------------------------------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------+ diff --git a/src/cython/doc/simplex_tree_sum.inc b/src/cython/doc/simplex_tree_sum.inc index fb0e54c1..086c69d5 100644 --- a/src/cython/doc/simplex_tree_sum.inc +++ b/src/cython/doc/simplex_tree_sum.inc @@ -1,14 +1,13 @@ -================================================================= =================================== =================================== -:Author: Clément Maria :Introduced in: GUDHI 2.0.0 :Copyright: GPL v3 -================================================================= =================================== =================================== +.. table:: + :widths: 30 50 20 -+----------------------------------------------------------------+------------------------------------------------------------------------+ -| .. figure:: | The simplex tree is an efficient and flexible data structure for | -| ../../doc/Simplex_tree/Simplex_tree_representation.png | representing general (filtered) simplicial complexes. | -| :alt: Simplex tree representation | | -| :figclass: align-center | The data structure is described in | -| | :cite:`boissonnatmariasimplextreealgorithmica` | -| Simplex tree representation | | -+----------------------------------------------------------------+------------------------------------------------------------------------+ -| :doc:`simplex_tree_user` | :doc:`simplex_tree_ref` | -+----------------------------------------------------------------+------------------------------------------------------------------------+ + +----------------------------------------------------------------+------------------------------------------------------------------------+-----------------------------+ + | .. figure:: | The simplex tree is an efficient and flexible data structure for | :Author: Clément Maria | + | ../../doc/Simplex_tree/Simplex_tree_representation.png | representing general (filtered) simplicial complexes. | | + | :alt: Simplex tree representation | | :Introduced in: GUDHI 2.0.0 | + | :figclass: align-center | The data structure is described in | | + | | :cite:`boissonnatmariasimplextreealgorithmica` | :Copyright: GPL v3 | + | | | | + +----------------------------------------------------------------+------------------------------------------------------------------------+-----------------------------+ + | * :doc:`simplex_tree_user` | * :doc:`simplex_tree_ref` | + +----------------------------------------------------------------+------------------------------------------------------------------------------------------------------+ diff --git a/src/cython/doc/tangential_complex_sum.inc b/src/cython/doc/tangential_complex_sum.inc index 72b4d7ba..0f03ffb3 100644 --- a/src/cython/doc/tangential_complex_sum.inc +++ b/src/cython/doc/tangential_complex_sum.inc @@ -1,15 +1,14 @@ -================================================================= =================================== =================================== -:Author: Clément Jamin :Introduced in: GUDHI 2.0.0 :Copyright: GPL v3 -:Requires: CGAL :math:`\geq` 4.8.0 Eigen3 -================================================================= =================================== =================================== +.. table:: + :widths: 30 50 20 -+----------------------------------------------------------------+------------------------------------------------------------------------+ -| .. figure:: | A Tangential Delaunay complex is a simplicial complex designed to | -| ../../doc/Tangential_complex/tc_examples.png | reconstruct a :math:`k`-dimensional manifold embedded in :math:`d`- | -| :figclass: align-center | dimensional Euclidean space. The input is a point sample coming from | -| | an unknown manifold. The running time depends only linearly on the | -| Tangential complex representation | extrinsic dimension :math:`d` and exponentially on the intrinsic | -| | dimension :math:`k`. | -+----------------------------------------------------------------+------------------------------------------------------------------------+ -| :doc:`tangential_complex_user` | :doc:`tangential_complex_ref` | -+----------------------------------------------------------------+------------------------------------------------------------------------+ + +----------------------------------------------------------------+------------------------------------------------------------------------+------------------------------------+ + | .. figure:: | A Tangential Delaunay complex is a simplicial complex designed to | :Author: Clément Jamin | + | ../../doc/Tangential_complex/tc_examples.png | reconstruct a :math:`k`-dimensional manifold embedded in :math:`d`- | | + | :figclass: align-center | dimensional Euclidean space. The input is a point sample coming from | :Introduced in: GUDHI 2.0.0 | + | | an unknown manifold. The running time depends only linearly on the | | + | | extrinsic dimension :math:`d` and exponentially on the intrinsic | :Copyright: GPL v3 | + | | dimension :math:`k`. | | + | | | :Requires: CGAL :math:`\geq` 4.8.0 | + +----------------------------------------------------------------+------------------------------------------------------------------------+------------------------------------+ + | * :doc:`tangential_complex_user` | * :doc:`tangential_complex_ref` | + +----------------------------------------------------------------+-------------------------------------------------------------------------------------------------------------+ diff --git a/src/cython/doc/witness_complex_sum.inc b/src/cython/doc/witness_complex_sum.inc index a8a126a0..49577745 100644 --- a/src/cython/doc/witness_complex_sum.inc +++ b/src/cython/doc/witness_complex_sum.inc @@ -1,19 +1,17 @@ -================================================================= =================================== =================================== -:Author: Siargey Kachanovich :Introduced in: GUDHI 2.0.0 :Copyright: GPL v3 -:Euclidean version requires: CGAL :math:`\geq` 4.6.0 Eigen3 -================================================================= =================================== =================================== +.. table:: + :widths: 30 50 20 -+-------------------------------------------------------------------+----------------------------------------------------------------------+ -| .. figure:: | Witness complex :math:`Wit(W,L)` is a simplicial complex defined on | -| ../../doc/Witness_complex/Witness_complex_representation.png | two sets of points in :math:`\mathbb{R}^D`. | -| :alt: Witness complex representation | | -| :figclass: align-center | The data structure is described in | -| | :cite:`boissonnatmariasimplextreealgorithmica`. | -| | | -| Witness complex representation | | -+-------------------------------------------------------------------+----------------------------------------------------------------------+ -| :doc:`witness_complex_user` | * :doc:`witness_complex_ref` | -| | * :doc:`strong_witness_complex_ref` | -| | * :doc:`euclidean_witness_complex_ref` | -| | * :doc:`euclidean_strong_witness_complex_ref` | -+-------------------------------------------------------------------+----------------------------------------------------------------------+ + +-------------------------------------------------------------------+----------------------------------------------------------------------+-----------------------------------------------------------------------------+ + | .. figure:: | Witness complex :math:`Wit(W,L)` is a simplicial complex defined on | :Author: Siargey Kachanovich | + | ../../doc/Witness_complex/Witness_complex_representation.png | two sets of points in :math:`\mathbb{R}^D`. | | + | :alt: Witness complex representation | | :Introduced in: GUDHI 2.0.0 | + | :figclass: align-center | The data structure is described in | | + | | :cite:`boissonnatmariasimplextreealgorithmica`. | :Copyright: GPL v3 | + | | | | + | | | :Requires: Eigen3 and CGAL :math:`\geq` 4.6.0 for Euclidean versions only | + +-------------------------------------------------------------------+----------------------------------------------------------------------+-----------------------------------------------------------------------------+ + | * :doc:`witness_complex_user` | * :doc:`witness_complex_ref` | + | | * :doc:`strong_witness_complex_ref` | + | | * :doc:`euclidean_witness_complex_ref` | + | | * :doc:`euclidean_strong_witness_complex_ref` | + +-------------------------------------------------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------+ diff --git a/src/cython/example/witness_complex_from_nearest_landmark_table.py b/src/cython/example/witness_complex_from_nearest_landmark_table.py index e6b295ee..1b79d9b2 100755 --- a/src/cython/example/witness_complex_from_nearest_landmark_table.py +++ b/src/cython/example/witness_complex_from_nearest_landmark_table.py @@ -30,14 +30,14 @@ __license__ = "GPL v3" print("#####################################################################") print("WitnessComplex creation from nearest landmark table") -nearest_landmark_table = [[[0, 0], [1, 1], [2, 2], [3, 3], [4, 4]], - [[1, 0], [2, 1], [3, 2], [4, 3], [0, 4]], - [[2, 0], [3, 1], [4, 2], [0, 3], [1, 4]], - [[3, 0], [4, 1], [0, 2], [1, 3], [2, 4]], - [[4, 0], [0, 1], [1, 2], [2, 3], [3, 4]]] +nearest_landmark_table = [[[0, 0.0], [1, 0.1], [2, 0.2], [3, 0.3], [4, 0.4]], + [[1, 0.0], [2, 0.1], [3, 0.2], [4, 0.3], [0, 0.4]], + [[2, 0.0], [3, 0.1], [4, 0.2], [0, 0.3], [1, 0.4]], + [[3, 0.0], [4, 0.1], [0, 0.2], [1, 0.3], [2, 0.4]], + [[4, 0.0], [0, 0.1], [1, 0.2], [2, 0.3], [3, 0.4]]] witness_complex = StrongWitnessComplex(nearest_landmark_table=nearest_landmark_table) -simplex_tree = witness_complex.create_simplex_tree(max_alpha_square=4.1) +simplex_tree = witness_complex.create_simplex_tree(max_alpha_square=0.41) message = "Number of simplices: " + repr(simplex_tree.num_simplices()) print(message) 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; |