diff options
Diffstat (limited to 'src/Bottleneck_distance')
23 files changed, 1709 insertions, 0 deletions
diff --git a/src/Bottleneck_distance/benchmark/CMakeLists.txt b/src/Bottleneck_distance/benchmark/CMakeLists.txt new file mode 100644 index 00000000..77cb013d --- /dev/null +++ b/src/Bottleneck_distance/benchmark/CMakeLists.txt @@ -0,0 +1,8 @@ +project(Bottleneck_distance_benchmark) + +if (NOT CGAL_VERSION VERSION_LESS 4.11.0) + add_executable ( bottleneck_chrono bottleneck_chrono.cpp ) + if (TBB_FOUND) + target_link_libraries(bottleneck_chrono ${TBB_LIBRARIES}) + endif(TBB_FOUND) +endif(NOT CGAL_VERSION VERSION_LESS 4.11.0) diff --git a/src/Bottleneck_distance/benchmark/bottleneck_chrono.cpp b/src/Bottleneck_distance/benchmark/bottleneck_chrono.cpp new file mode 100644 index 00000000..576d510b --- /dev/null +++ b/src/Bottleneck_distance/benchmark/bottleneck_chrono.cpp @@ -0,0 +1,50 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author: Francois Godi + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <gudhi/Bottleneck.h> +#include <chrono> +#include <fstream> +#include <random> + +using namespace Gudhi::persistence_diagram; + + +double upper_bound = 400.; // any real > 0 + +int main() { + std::ofstream result_file; + result_file.open("results.csv", std::ios::out); + + for (int n = 1000; n <= 10000; n += 1000) { + std::uniform_real_distribution<double> unif1(0., upper_bound); + std::uniform_real_distribution<double> unif2(upper_bound / 1000., upper_bound / 100.); + std::default_random_engine re; + std::vector< std::pair<double, double> > v1, v2; + for (int i = 0; i < n; i++) { + double a = unif1(re); + double b = unif1(re); + double x = unif2(re); + double y = unif2(re); + v1.emplace_back(std::min(a, b), std::max(a, b)); + v2.emplace_back(std::min(a, b) + std::min(x, y), std::max(a, b) + std::max(x, y)); + if (i % 5 == 0) + v1.emplace_back(std::min(a, b), std::min(a, b) + x); + if (i % 3 == 0) + v2.emplace_back(std::max(a, b), std::max(a, b) + y); + } + std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now(); + double b = bottleneck_distance(v1, v2); + std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now(); + typedef std::chrono::duration<int, std::milli> millisecs_t; + millisecs_t duration(std::chrono::duration_cast<millisecs_t>(end - start)); + result_file << n << ";" << duration.count() << ";" << b << std::endl; + } + result_file.close(); +} diff --git a/src/Bottleneck_distance/concept/Persistence_diagram.h b/src/Bottleneck_distance/concept/Persistence_diagram.h new file mode 100644 index 00000000..8c8761cb --- /dev/null +++ b/src/Bottleneck_distance/concept/Persistence_diagram.h @@ -0,0 +1,38 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author: François Godi + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef CONCEPT_BOTTLENECK_DISTANCE_PERSISTENCE_DIAGRAM_H_ +#define CONCEPT_BOTTLENECK_DISTANCE_PERSISTENCE_DIAGRAM_H_ + +namespace Gudhi { + +namespace persistence_diagram { + +/** \brief Concept of point in a persistence diagram. std::get<0>(point) must return the birth of the corresponding component and std::get<1>(point) its death. + * Both should be convertible to `double`. + * A valid implementation of this concept is std::pair<double,double>. + * Death should be larger than birth, death can be std::numeric_limits<double>::infinity() for components which stay alive. + * + * \ingroup bottleneck_distance + */ +struct DiagramPoint{}; + +/** \brief Concept of persistence diagram. It is a range of `DiagramPoint`. + * std::begin(diagram) and std::end(diagram) must return corresponding iterators. + * + * \ingroup bottleneck_distance + */ +struct PersistenceDiagram{}; + +} // namespace persistence_diagram + +} // namespace Gudhi + +#endif // CONCEPT_BOTTLENECK_DISTANCE_PERSISTENCE_DIAGRAM_H_ diff --git a/src/Bottleneck_distance/doc/COPYRIGHT b/src/Bottleneck_distance/doc/COPYRIGHT new file mode 100644 index 00000000..61f17f6d --- /dev/null +++ b/src/Bottleneck_distance/doc/COPYRIGHT @@ -0,0 +1,12 @@ +The files of this directory are part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. +See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + +Author(s): Vincent Rouvreau + +Copyright (C) 2015 Inria + +This gives everyone the freedoms to use openFrameworks in any context: +commercial or non-commercial, public or private, open or closed source. + +You should have received a copy of the MIT License along with this program. +If not, see https://opensource.org/licenses/MIT.
\ No newline at end of file diff --git a/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h b/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h new file mode 100644 index 00000000..bbc952e1 --- /dev/null +++ b/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h @@ -0,0 +1,81 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author: François Godi + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + * - 2019/06 Vincent Rouvreau : Fix #11 - Distance computation shall be better documented. + */ + +#ifndef DOC_BOTTLENECK_DISTANCE_INTRO_BOTTLENECK_DISTANCE_H_ +#define DOC_BOTTLENECK_DISTANCE_INTRO_BOTTLENECK_DISTANCE_H_ + +// needs namespace for Doxygen to link on classes +namespace Gudhi { +// needs namespace for Doxygen to link on classes +namespace persistence_diagram { + +/** \defgroup bottleneck_distance Bottleneck distance + * + * \author François Godi + * @{ + * + * \section bottleneckdefinition Definition + * + * The bottleneck distance measures the similarity between two persistence diagrams. It is the shortest distance b for + * which there exists a perfect matching between the points of the two diagrams (completed with all the points on the + * diagonal in order to ignore cardinality mismatchs) such that any couple of matched points are at distance at most b, + * where the distance between points is the sup norm in \f$\mathbb{R}^2\f$ (not the Euclidean distance). + * + * \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. + * + * \section bottleneckdistanceprecision Distance computation + * + * The following example explains how the distance is computed: + * + * \code{.cpp} +#include <gudhi/Bottleneck.h> + +#include <iostream> +#include <vector> +#include <utility> // for pair + +int main() { + std::vector< std::pair<double, double> > diag1, diag2; + diag1.emplace_back(0., 0.); + diag2.emplace_back(0., 13.); + + double b = Gudhi::persistence_diagram::bottleneck_distance(diag1, diag2); + std::cout << "Bottleneck distance = " << b << std::endl; +} + * \endcode + * + * \code Bottleneck distance = 6.5 + * \endcode + * + * \image html bottleneck_distance_example.png The point (0, 13) is at distance 6.5 from the diagonal and more specifically from the point (6.5, 6.5) + * + * \section bottleneckbasicexample Basic example + * + * This other example computes the bottleneck distance from 2 persistence diagrams: + * \include Bottleneck_distance/bottleneck_basic_example.cpp + * + * \code + Bottleneck distance = 0.75 + Approx bottleneck distance = 0.808176 + * \endcode + + */ +/** @} */ // end defgroup bottleneck_distance + +} // namespace persistence_diagram + +} // namespace Gudhi + +#endif // DOC_BOTTLENECK_DISTANCE_INTRO_BOTTLENECK_DISTANCE_H_ diff --git a/src/Bottleneck_distance/doc/bottleneck_distance_example.ipe b/src/Bottleneck_distance/doc/bottleneck_distance_example.ipe new file mode 100644 index 00000000..2033ea56 --- /dev/null +++ b/src/Bottleneck_distance/doc/bottleneck_distance_example.ipe @@ -0,0 +1,287 @@ +<?xml version="1.0"?> +<!DOCTYPE ipe SYSTEM "ipe.dtd"> +<ipe version="70206" creator="Ipe 7.2.7"> +<info created="D:20190606115351" modified="D:20190607172542"/> +<ipestyle name="basic"> +<symbol name="arrow/arc(spx)"> +<path stroke="sym-stroke" fill="sym-stroke" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/farc(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/ptarc(spx)"> +<path stroke="sym-stroke" fill="sym-stroke" pen="sym-pen"> +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/fptarc(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="mark/circle(sx)" transformations="translations"> +<path fill="sym-stroke"> +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e +</path> +</symbol> +<symbol name="mark/disk(sx)" transformations="translations"> +<path fill="sym-stroke"> +0.6 0 0 0.6 0 0 e +</path> +</symbol> +<symbol name="mark/fdisk(sfx)" transformations="translations"> +<group> +<path fill="sym-fill"> +0.5 0 0 0.5 0 0 e +</path> +<path fill="sym-stroke" fillrule="eofill"> +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e +</path> +</group> +</symbol> +<symbol name="mark/box(sx)" transformations="translations"> +<path fill="sym-stroke" fillrule="eofill"> +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +-0.4 -0.4 m +0.4 -0.4 l +0.4 0.4 l +-0.4 0.4 l +h +</path> +</symbol> +<symbol name="mark/square(sx)" transformations="translations"> +<path fill="sym-stroke"> +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +</path> +</symbol> +<symbol name="mark/fsquare(sfx)" transformations="translations"> +<group> +<path fill="sym-fill"> +-0.5 -0.5 m +0.5 -0.5 l +0.5 0.5 l +-0.5 0.5 l +h +</path> +<path fill="sym-stroke" fillrule="eofill"> +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +-0.4 -0.4 m +0.4 -0.4 l +0.4 0.4 l +-0.4 0.4 l +h +</path> +</group> +</symbol> +<symbol name="mark/cross(sx)" transformations="translations"> +<group> +<path fill="sym-stroke"> +-0.43 -0.57 m +0.57 0.43 l +0.43 0.57 l +-0.57 -0.43 l +h +</path> +<path fill="sym-stroke"> +-0.43 0.57 m +0.57 -0.43 l +0.43 -0.57 l +-0.57 0.43 l +h +</path> +</group> +</symbol> +<symbol name="arrow/fnormal(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/pointed(spx)"> +<path stroke="sym-stroke" fill="sym-stroke" pen="sym-pen"> +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/fpointed(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/linear(spx)"> +<path stroke="sym-stroke" pen="sym-pen"> +-1 0.333 m +0 0 l +-1 -0.333 l +</path> +</symbol> +<symbol name="arrow/fdouble(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/double(spx)"> +<path stroke="sym-stroke" fill="sym-stroke" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h +</path> +</symbol> +<pen name="heavier" value="0.8"/> +<pen name="fat" value="1.2"/> +<pen name="ultrafat" value="2"/> +<symbolsize name="large" value="5"/> +<symbolsize name="small" value="2"/> +<symbolsize name="tiny" value="1.1"/> +<arrowsize name="large" value="10"/> +<arrowsize name="small" value="5"/> +<arrowsize name="tiny" value="3"/> +<color name="red" value="1 0 0"/> +<color name="green" value="0 1 0"/> +<color name="blue" value="0 0 1"/> +<color name="yellow" value="1 1 0"/> +<color name="orange" value="1 0.647 0"/> +<color name="gold" value="1 0.843 0"/> +<color name="purple" value="0.627 0.125 0.941"/> +<color name="gray" value="0.745"/> +<color name="brown" value="0.647 0.165 0.165"/> +<color name="navy" value="0 0 0.502"/> +<color name="pink" value="1 0.753 0.796"/> +<color name="seagreen" value="0.18 0.545 0.341"/> +<color name="turquoise" value="0.251 0.878 0.816"/> +<color name="violet" value="0.933 0.51 0.933"/> +<color name="darkblue" value="0 0 0.545"/> +<color name="darkcyan" value="0 0.545 0.545"/> +<color name="darkgray" value="0.663"/> +<color name="darkgreen" value="0 0.392 0"/> +<color name="darkmagenta" value="0.545 0 0.545"/> +<color name="darkorange" value="1 0.549 0"/> +<color name="darkred" value="0.545 0 0"/> +<color name="lightblue" value="0.678 0.847 0.902"/> +<color name="lightcyan" value="0.878 1 1"/> +<color name="lightgray" value="0.827"/> +<color name="lightgreen" value="0.565 0.933 0.565"/> +<color name="lightyellow" value="1 1 0.878"/> +<dashstyle name="dashed" value="[4] 0"/> +<dashstyle name="dotted" value="[1 3] 0"/> +<dashstyle name="dash dotted" value="[4 2 1 2] 0"/> +<dashstyle name="dash dot dotted" value="[4 2 1 2 1 2] 0"/> +<textsize name="large" value="\large"/> +<textsize name="small" value="\small"/> +<textsize name="tiny" value="\tiny"/> +<textsize name="Large" value="\Large"/> +<textsize name="LARGE" value="\LARGE"/> +<textsize name="huge" value="\huge"/> +<textsize name="Huge" value="\Huge"/> +<textsize name="footnote" value="\footnotesize"/> +<textstyle name="center" begin="\begin{center}" end="\end{center}"/> +<textstyle name="itemize" begin="\begin{itemize}" end="\end{itemize}"/> +<textstyle name="item" begin="\begin{itemize}\item{}" end="\end{itemize}"/> +<gridsize name="4 pts" value="4"/> +<gridsize name="8 pts (~3 mm)" value="8"/> +<gridsize name="16 pts (~6 mm)" value="16"/> +<gridsize name="32 pts (~12 mm)" value="32"/> +<gridsize name="10 pts (~3.5 mm)" value="10"/> +<gridsize name="20 pts (~7 mm)" value="20"/> +<gridsize name="14 pts (~5 mm)" value="14"/> +<gridsize name="28 pts (~10 mm)" value="28"/> +<gridsize name="56 pts (~20 mm)" value="56"/> +<anglesize name="90 deg" value="90"/> +<anglesize name="60 deg" value="60"/> +<anglesize name="45 deg" value="45"/> +<anglesize name="30 deg" value="30"/> +<anglesize name="22.5 deg" value="22.5"/> +<opacity name="10%" value="0.1"/> +<opacity name="30%" value="0.3"/> +<opacity name="50%" value="0.5"/> +<opacity name="75%" value="0.75"/> +<tiling name="falling" angle="-60" step="4" width="1"/> +<tiling name="rising" angle="30" step="4" width="1"/> +</ipestyle> +<page> +<layer name="alpha"/> +<view layers="alpha" active="alpha"/> +<path layer="alpha" stroke="black" arrow="normal/normal"> +160 256 m +288 256 l +</path> +<path stroke="black" arrow="normal/normal"> +160 256 m +160 384 l +</path> +<path stroke="navy"> +160 256 m +264 360 l +264 360 l +</path> +<use name="mark/fdisk(sfx)" pos="160 256" size="normal" stroke="red" fill="red"/> +<use name="mark/fdisk(sfx)" pos="160 360" size="normal" stroke="navy" fill="navy"/> +<path stroke="red"> +200 296 m +224 320 l +</path> +<text matrix="1 0 0 1 -32 8" transformations="translations" pos="164 248" stroke="red" type="label" width="23.8" height="7.473" depth="2.49" valign="baseline">(0, 0)</text> +<text matrix="1 0 0 1 -4 4" transformations="translations" pos="132 356" stroke="navy" type="label" width="28.781" height="7.473" depth="2.49" valign="baseline">(0, 13)</text> +<text transformations="translations" pos="216 300" stroke="black" type="label" width="39.297" height="7.473" depth="2.49" valign="baseline">(6.5, 6.5)</text> +<path stroke="black" arrow="normal/normal" rarrow="normal/normal"> +160.433 359.995 m +213.999 309.986 l +</path> +<use name="mark/cross(sx)" pos="214.209 309.567" size="normal" stroke="black"/> +<text matrix="1 0 0 1 6.00021 -25.5247" transformations="translations" pos="174.375 385.835" stroke="black" type="minipage" width="73.2406" height="11.924" depth="6.95" valign="top">Bottleneck +distance is 6.5</text> +</page> +</ipe> diff --git a/src/Bottleneck_distance/doc/bottleneck_distance_example.png b/src/Bottleneck_distance/doc/bottleneck_distance_example.png Binary files differnew file mode 100644 index 00000000..1d3b91aa --- /dev/null +++ b/src/Bottleneck_distance/doc/bottleneck_distance_example.png diff --git a/src/Bottleneck_distance/doc/perturb_pd.png b/src/Bottleneck_distance/doc/perturb_pd.png Binary files differnew file mode 100644 index 00000000..be638de0 --- /dev/null +++ b/src/Bottleneck_distance/doc/perturb_pd.png diff --git a/src/Bottleneck_distance/example/CMakeLists.txt b/src/Bottleneck_distance/example/CMakeLists.txt new file mode 100644 index 00000000..3d65963a --- /dev/null +++ b/src/Bottleneck_distance/example/CMakeLists.txt @@ -0,0 +1,25 @@ +project(Bottleneck_distance_examples) + +if (NOT CGAL_VERSION VERSION_LESS 4.11.0) + add_executable (bottleneck_basic_example bottleneck_basic_example.cpp) + + if (TBB_FOUND) + target_link_libraries(bottleneck_basic_example ${TBB_LIBRARIES}) + endif(TBB_FOUND) + + add_test(NAME Bottleneck_distance_example_basic COMMAND $<TARGET_FILE:bottleneck_basic_example>) + +endif (NOT CGAL_VERSION VERSION_LESS 4.11.0) + +if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + add_executable (alpha_rips_persistence_bottleneck_distance alpha_rips_persistence_bottleneck_distance.cpp) + target_link_libraries(alpha_rips_persistence_bottleneck_distance ${Boost_PROGRAM_OPTIONS_LIBRARY}) + + if (TBB_FOUND) + target_link_libraries(alpha_rips_persistence_bottleneck_distance ${TBB_LIBRARIES}) + endif(TBB_FOUND) + + add_test(NAME Bottleneck_distance_example_alpha_rips_persistence_bottleneck + COMMAND $<TARGET_FILE:alpha_rips_persistence_bottleneck_distance> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.15" "-m" "0.12" "-d" "3" "-p" "3") +endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) diff --git a/src/Bottleneck_distance/example/README b/src/Bottleneck_distance/example/README new file mode 100644 index 00000000..01bcd74a --- /dev/null +++ b/src/Bottleneck_distance/example/README @@ -0,0 +1,19 @@ +# Bottleneck_distance #
+
+## `alpha_rips_persistence_bottleneck_distance` ##
+This program computes the persistent homology with coefficient field Z/pZ of a Rips complex defined on a set of input points. The output diagram contains one bar per line, written with the convention:
+
+`p dim birth death`
+
+where `dim` is the dimension of the homological feature, `birth` and `death` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients.
+
+Usage:
+`alpha_rips_persistence_bottleneck_distance [options] <OFF input file>`
+
+Allowed options:
+
+* `-h [ --help ]` Produce help message
+* `-r [ --max-edge-length ]` (default = inf) Maximal length of an edge for the Rips complex construction.`
+* `-d [ --cpx-dimension ]` (default = 1) Maximal dimension of the Rips complex we want to compute.`
+* `-p [ --field-charac ]` (default = 11) Characteristic p of the coefficient field Z/pZ for computing homology.
+* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals.
diff --git a/src/Bottleneck_distance/example/alpha_rips_persistence_bottleneck_distance.cpp b/src/Bottleneck_distance/example/alpha_rips_persistence_bottleneck_distance.cpp new file mode 100644 index 00000000..6c0dc9bf --- /dev/null +++ b/src/Bottleneck_distance/example/alpha_rips_persistence_bottleneck_distance.cpp @@ -0,0 +1,178 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2017 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <gudhi/Alpha_complex.h> +#include <gudhi/Rips_complex.h> +#include <gudhi/distance_functions.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Points_off_io.h> +#include <gudhi/Bottleneck.h> + +#include <CGAL/Epick_d.h> + +#include <boost/program_options.hpp> + +#include <string> +#include <vector> +#include <limits> // infinity +#include <utility> // for pair +#include <algorithm> // for transform + + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp >; +using Kernel = CGAL::Epick_d< CGAL::Dynamic_dimension_tag >; +using Point_d = Kernel::Point_d; +using Points_off_reader = Gudhi::Points_off_reader<Point_d>; + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , Filtration_value & threshold + , int & dim_max + , int & p + , Filtration_value & min_persistence); + +static inline std::pair<double, double> compute_root_square(std::pair<double, double> input) { + return std::make_pair(std::sqrt(input.first), std::sqrt(input.second)); +} + +int main(int argc, char * argv[]) { + std::string off_file_points; + Filtration_value threshold; + int dim_max; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, threshold, dim_max, p, min_persistence); + + Points_off_reader off_reader(off_file_points); + + // -------------------------------------------- + // Rips persistence + // -------------------------------------------- + Rips_complex rips_complex(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); + + // Construct the Rips complex in a Simplex Tree + Simplex_tree rips_stree; + + rips_complex.create_complex(rips_stree, dim_max); + std::cout << "The Rips complex contains " << rips_stree.num_simplices() << " simplices and has dimension " + << rips_stree.dimension() << " \n"; + + // Sort the simplices in the order of the filtration + rips_stree.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology rips_pcoh(rips_stree); + // initializes the coefficient field for homology + rips_pcoh.init_coefficients(p); + rips_pcoh.compute_persistent_cohomology(min_persistence); + + // rips_pcoh.output_diagram(); + + // -------------------------------------------- + // Alpha persistence + // -------------------------------------------- + Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex(off_reader.get_point_cloud()); + + Simplex_tree alpha_stree; + alpha_complex.create_complex(alpha_stree, threshold * threshold); + std::cout << "The Alpha complex contains " << alpha_stree.num_simplices() << " simplices and has dimension " + << alpha_stree.dimension() << " \n"; + + // Sort the simplices in the order of the filtration + alpha_stree.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology alpha_pcoh(alpha_stree); + // initializes the coefficient field for homology + alpha_pcoh.init_coefficients(p); + alpha_pcoh.compute_persistent_cohomology(min_persistence * min_persistence); + + // alpha_pcoh.output_diagram(); + + // -------------------------------------------- + // Bottleneck distance between both persistence + // -------------------------------------------- + double max_b_distance {}; + for (int dim = 0; dim < dim_max; dim ++) { + std::vector< std::pair< Filtration_value , Filtration_value > > rips_intervals; + std::vector< std::pair< Filtration_value , Filtration_value > > alpha_intervals; + rips_intervals = rips_pcoh.intervals_in_dimension(dim); + alpha_intervals = alpha_pcoh.intervals_in_dimension(dim); + std::transform(alpha_intervals.begin(), alpha_intervals.end(), alpha_intervals.begin(), compute_root_square); + + double bottleneck_distance = Gudhi::persistence_diagram::bottleneck_distance(rips_intervals, alpha_intervals); + std::cout << "In dimension " << dim << ", bottleneck distance = " << bottleneck_distance << std::endl; + if (bottleneck_distance > max_b_distance) + max_b_distance = bottleneck_distance; + } + std::cout << "================================================================================" << std::endl; + std::cout << "Bottleneck distance is " << max_b_distance << std::endl; + + return 0; +} + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , Filtration_value & threshold + , int & dim_max + , int & p + , Filtration_value & min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options() + ("input-file", po::value<std::string>(&off_file_points), + "Name of an OFF file containing a point set.\n"); + + po::options_description visible("Allowed options", 100); + visible.add_options() + ("help,h", "produce help message") + ("max-edge-length,r", + po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()), + "Maximal length of an edge for the Rips complex construction.") + ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), + "Maximal dimension of the Rips complex we want to compute.") + ("field-charac,p", po::value<int>(&p)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.") + ("min-persistence,m", po::value<Filtration_value>(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals"); + + po::positional_options_description pos; + pos.add("input-file", 1); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map vm; + po::store(po::command_line_parser(argc, argv). + options(all).positional(pos).run(), vm); + po::notify(vm); + + if (vm.count("help") || !vm.count("input-file")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a Rips complex defined on a set of input points.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + exit(-1); + } +} diff --git a/src/Bottleneck_distance/example/bottleneck_basic_example.cpp b/src/Bottleneck_distance/example/bottleneck_basic_example.cpp new file mode 100644 index 00000000..61778a55 --- /dev/null +++ b/src/Bottleneck_distance/example/bottleneck_basic_example.cpp @@ -0,0 +1,28 @@ +#include <gudhi/Bottleneck.h> + +#include <iostream> +#include <vector> +#include <utility> // for pair +#include <limits> // for numeric_limits + +int main() { + std::vector< std::pair<double, double> > v1, v2; + + v1.emplace_back(2.7, 3.7); + v1.emplace_back(9.6, 14.); + v1.emplace_back(34.2, 34.974); + v1.emplace_back(3., std::numeric_limits<double>::infinity()); + + v2.emplace_back(2.8, 4.45); + v2.emplace_back(9.5, 14.1); + v2.emplace_back(3.2, std::numeric_limits<double>::infinity()); + + + double b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2); + + std::cout << "Bottleneck distance = " << b << std::endl; + + b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2, 0.1); + + std::cout << "Approx bottleneck distance = " << b << std::endl; +} diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h new file mode 100644 index 00000000..82ba9f68 --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h @@ -0,0 +1,124 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author: Francois Godi + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - 2019/06 Vincent Rouvreau : Fix doxygen warning. + * - 2019/08 Vincent Rouvreau: Fix issue #10 for CGAL + * - YYYY/MM Author: Description of the modification + */ + +#ifndef BOTTLENECK_H_ +#define BOTTLENECK_H_ + +#include <gudhi/Graph_matching.h> + +#include <CGAL/version.h> // for CGAL_VERSION_NR + +#include <vector> +#include <algorithm> // for max +#include <limits> // for numeric_limits + +#include <cmath> +#include <cfloat> // FLT_EVAL_METHOD + +// Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10 +#if CGAL_VERSION_NR < 1041101000 +# error Alpha_complex_3d is only available for CGAL >= 4.11 +#endif + +namespace Gudhi { + +namespace persistence_diagram { + +inline double bottleneck_distance_approx(Persistence_graph& g, double e) { + double b_lower_bound = 0.; + double b_upper_bound = g.diameter_bound(); + const double alpha = std::pow(g.size(), 1. / 5.); + Graph_matching m(g); + Graph_matching biggest_unperfect(g); + while (b_upper_bound - b_lower_bound > 2 * e) { + double step = b_lower_bound + (b_upper_bound - b_lower_bound) / alpha; +#if !defined FLT_EVAL_METHOD || FLT_EVAL_METHOD < 0 || FLT_EVAL_METHOD > 1 + // On platforms where double computation is done with excess precision, + // we force it to its true precision so the following test is reliable. + volatile double drop_excess_precision = step; + step = drop_excess_precision; + // Alternative: step = CGAL::IA_force_to_double(step); +#endif + if (step <= b_lower_bound || step >= b_upper_bound) // Avoid precision problem + break; + m.set_r(step); + while (m.multi_augment()) {} // compute a maximum matching (in the graph corresponding to the current r) + if (m.perfect()) { + m = biggest_unperfect; + b_upper_bound = step; + } else { + biggest_unperfect = m; + b_lower_bound = step; + } + } + return (b_lower_bound + b_upper_bound) / 2.; +} + +inline double bottleneck_distance_exact(Persistence_graph& g) { + std::vector<double> sd = g.sorted_distances(); + long lower_bound_i = 0; + long upper_bound_i = sd.size() - 1; + const double alpha = std::pow(g.size(), 1. / 5.); + Graph_matching m(g); + Graph_matching biggest_unperfect(g); + while (lower_bound_i != upper_bound_i) { + long step = lower_bound_i + static_cast<long> ((upper_bound_i - lower_bound_i - 1) / alpha); + m.set_r(sd.at(step)); + while (m.multi_augment()) {} // compute a maximum matching (in the graph corresponding to the current r) + if (m.perfect()) { + m = biggest_unperfect; + upper_bound_i = step; + } else { + biggest_unperfect = m; + lower_bound_i = step + 1; + } + } + return sd.at(lower_bound_i); +} + +/** \brief Function to compute the Bottleneck distance between two persistence diagrams. + * + * \tparam Persistence_diagram1,Persistence_diagram2 + * models of the concept `PersistenceDiagram`. + * + * \param[in] diag1 The first persistence diagram. + * \param[in] diag2 The second persistence diagram. + * + * \param[in] e + * \parblock + * If `e` is 0, this uses an expensive algorithm to compute the exact distance. + * + * If `e` is not 0, it asks for an additive `e`-approximation, and currently + * also allows a small multiplicative error (the last 2 or 3 bits of the + * mantissa may be wrong). This version of the algorithm takes advantage of the + * limited precision of `double` and is usually a lot faster to compute, + * whatever the value of `e`. + * + * Thus, by default, `e` is the smallest positive double. + * \endparblock + * + * \ingroup bottleneck_distance + */ +template<typename Persistence_diagram1, typename Persistence_diagram2> +double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, + double e = (std::numeric_limits<double>::min)()) { + Persistence_graph g(diag1, diag2, e); + if (g.bottleneck_alive() == std::numeric_limits<double>::infinity()) + return std::numeric_limits<double>::infinity(); + return (std::max)(g.bottleneck_alive(), e == 0. ? bottleneck_distance_exact(g) : bottleneck_distance_approx(g, e)); +} + +} // namespace persistence_diagram + +} // namespace Gudhi + +#endif // BOTTLENECK_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Graph_matching.h b/src/Bottleneck_distance/include/gudhi/Graph_matching.h new file mode 100644 index 00000000..7b7623ce --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Graph_matching.h @@ -0,0 +1,162 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author: Francois Godi + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef GRAPH_MATCHING_H_ +#define GRAPH_MATCHING_H_ + +#include <gudhi/Neighbors_finder.h> + +#include <vector> +#include <unordered_set> +#include <algorithm> + +namespace Gudhi { + +namespace persistence_diagram { + +/** \internal \brief Structure representing a graph matching. The graph is a Persistence_diagrams_graph. + * + * \ingroup bottleneck_distance + */ +class Graph_matching { + public: + /** \internal \brief Constructor constructing an empty matching. */ + explicit Graph_matching(Persistence_graph &g); + /** \internal \brief Is the matching perfect ? */ + bool perfect() const; + /** \internal \brief Augments the matching with a maximal set of edge-disjoint shortest augmenting paths. */ + bool multi_augment(); + /** \internal \brief Sets the maximum length of the edges allowed to be added in the matching, 0 initially. */ + void set_r(double r); + + private: + Persistence_graph* gp; + double r; + /** \internal \brief Given a point from V, provides its matched point in U, null_point_index() if there isn't. */ + std::vector<int> v_to_u; + /** \internal \brief All the unmatched points in U. */ + std::unordered_set<int> unmatched_in_u; + + /** \internal \brief Provides a Layered_neighbors_finder dividing the graph in layers. Basically a BFS. */ + Layered_neighbors_finder layering() const; + /** \internal \brief Augments the matching with a simple path no longer than max_depth. Basically a DFS. */ + bool augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth); + /** \internal \brief Update the matching with the simple augmenting path given as parameter. */ + void update(std::vector<int> & path); +}; + +inline Graph_matching::Graph_matching(Persistence_graph& g) + : gp(&g), r(0.), v_to_u(g.size(), null_point_index()), unmatched_in_u(g.size()) { + for (int u_point_index = 0; u_point_index < g.size(); ++u_point_index) + unmatched_in_u.insert(u_point_index); +} + +inline bool Graph_matching::perfect() const { + return unmatched_in_u.empty(); +} + +inline bool Graph_matching::multi_augment() { + if (perfect()) + return false; + Layered_neighbors_finder layered_nf(layering()); + int max_depth = layered_nf.vlayers_number()*2 - 1; + double rn = sqrt(gp->size()); + // verification of a necessary criterion in order to shortcut if possible + if (max_depth < 0 || (unmatched_in_u.size() > rn && max_depth >= rn)) + return false; + bool successful = false; + std::vector<int> tries(unmatched_in_u.cbegin(), unmatched_in_u.cend()); + for (auto it = tries.cbegin(); it != tries.cend(); it++) + // 'augment' has side-effects which have to be always executed, don't change order + successful = augment(layered_nf, *it, max_depth) || successful; + return successful; +} + +inline void Graph_matching::set_r(double r) { + this->r = r; +} + +inline bool Graph_matching::augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth) { + // V vertices have at most one successor, thus when we backtrack from U we can directly pop_back 2 vertices. + std::vector<int> path; + path.emplace_back(u_start_index); + do { + if (static_cast<int> (path.size()) > max_depth) { + path.pop_back(); + path.pop_back(); + } + if (path.empty()) + return false; + path.emplace_back(layered_nf.pull_near(path.back(), static_cast<int> (path.size()) / 2)); + while (path.back() == null_point_index()) { + path.pop_back(); + path.pop_back(); + if (path.empty()) + return false; + path.pop_back(); + path.emplace_back(layered_nf.pull_near(path.back(), path.size() / 2)); + } + path.emplace_back(v_to_u.at(path.back())); + } while (path.back() != null_point_index()); + // if v_to_u.at(path.back()) has no successor, path.back() is an exposed vertex + path.pop_back(); + update(path); + return true; +} + +inline Layered_neighbors_finder Graph_matching::layering() const { + std::vector<int> u_vertices(unmatched_in_u.cbegin(), unmatched_in_u.cend()); + std::vector<int> v_vertices; + Neighbors_finder nf(*gp, r); + for (int v_point_index = 0; v_point_index < gp->size(); ++v_point_index) + nf.add(v_point_index); + Layered_neighbors_finder layered_nf(*gp, r); + for (int layer = 0; !u_vertices.empty(); layer++) { + // one layer is one step in the BFS + for (auto it1 = u_vertices.cbegin(); it1 != u_vertices.cend(); ++it1) { + std::vector<int> u_succ(nf.pull_all_near(*it1)); + for (auto it2 = u_succ.begin(); it2 != u_succ.end(); ++it2) { + layered_nf.add(*it2, layer); + v_vertices.emplace_back(*it2); + } + } + // When the above for finishes, we have progress of one half-step (from U to V) in the BFS + u_vertices.clear(); + bool end = false; + for (auto it = v_vertices.cbegin(); it != v_vertices.cend(); it++) + if (v_to_u.at(*it) == null_point_index()) + // we stop when a nearest exposed V vertex (from U exposed vertices) has been found + end = true; + else + u_vertices.emplace_back(v_to_u.at(*it)); + // When the above for finishes, we have progress of one half-step (from V to U) in the BFS + if (end) + return layered_nf; + v_vertices.clear(); + } + return layered_nf; +} + +inline void Graph_matching::update(std::vector<int>& path) { + // Must return 1. + unmatched_in_u.erase(path.front()); + for (auto it = path.cbegin(); it != path.cend(); ++it) { + // Be careful, the iterator is incremented twice each time + int tmp = *it; + v_to_u[*(++it)] = tmp; + } +} + + +} // namespace persistence_diagram + +} // namespace Gudhi + +#endif // GRAPH_MATCHING_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Internal_point.h b/src/Bottleneck_distance/include/gudhi/Internal_point.h new file mode 100644 index 00000000..f829943e --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Internal_point.h @@ -0,0 +1,79 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author: Francois Godi + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef INTERNAL_POINT_H_ +#define INTERNAL_POINT_H_ + +namespace Gudhi { + +namespace persistence_diagram { + +/** \internal \brief Returns the used index for encoding none of the points */ +int null_point_index(); + +/** \internal \typedef \brief Internal_point is the internal points representation, indexes used outside. */ +struct Internal_point { + double vec[2]; + int point_index; + + Internal_point() { } + + Internal_point(double x, double y, int p_i) { + vec[0] = x; + vec[1] = y; + point_index = p_i; + } + + double x() const { + return vec[ 0 ]; + } + + double y() const { + return vec[ 1 ]; + } + + double& x() { + return vec[ 0 ]; + } + + double& y() { + return vec[ 1 ]; + } + + bool operator==(const Internal_point& p) const { + return point_index == p.point_index; + } + + bool operator!=(const Internal_point& p) const { + return !(*this == p); + } +}; + +inline int null_point_index() { + return -1; +} + +struct Construct_coord_iterator { + typedef const double* result_type; + + const double* operator()(const Internal_point& p) const { + return p.vec; + } + + const double* operator()(const Internal_point& p, int) const { + return p.vec + 2; + } +}; + +} // namespace persistence_diagram + +} // namespace Gudhi + +#endif // INTERNAL_POINT_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h new file mode 100644 index 00000000..c65e6082 --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h @@ -0,0 +1,182 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author: Francois Godi + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef NEIGHBORS_FINDER_H_ +#define NEIGHBORS_FINDER_H_ + +// Inclusion order is important for CGAL patch +#include <CGAL/Kd_tree.h> +#include <CGAL/Search_traits.h> + +#include <gudhi/Persistence_graph.h> +#include <gudhi/Internal_point.h> + +#include <unordered_set> +#include <vector> +#include <algorithm> // for std::max +#include <cmath> // for std::abs + +namespace Gudhi { + +namespace persistence_diagram { + +/** \internal \brief Variant of CGAL::Fuzzy_iso_box to ensure that the box ic closed. It isn't clear how necessary that is. + */ +struct Square_query { + typedef CGAL::Dimension_tag<2> D; + 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; + } + bool inner_range_intersects(CGAL::Kd_tree_rectangle<FT, D> const&r) const { + return + r.max_coord(0) >= c.x() - size && + r.min_coord(0) <= c.x() + size && + r.max_coord(1) >= c.y() - size && + r.min_coord(1) <= c.y() + size; + } + bool outer_range_contains(CGAL::Kd_tree_rectangle<FT, D> const&r) const { + return + r.min_coord(0) >= c.x() - size && + r.max_coord(0) <= c.x() + size && + r.min_coord(1) >= c.y() - size && + r.max_coord(1) <= c.y() + size; + } + Point_d c; + FT size; +}; + +/** \internal \brief data structure used to find any point (including projections) in V near to a query point from U + * (which can be a projection). + * + * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically + * removed. + * + * \ingroup bottleneck_distance + */ +class Neighbors_finder { + typedef CGAL::Dimension_tag<2> D; + typedef CGAL::Search_traits<double, Internal_point, const double*, Construct_coord_iterator, D> Traits; + typedef CGAL::Kd_tree<Traits> Kd_tree; + + public: + /** \internal \brief Constructor taking the near distance definition as parameter. */ + Neighbors_finder(const Persistence_graph& g, double r); + /** \internal \brief A point added will be possibly pulled. */ + void add(int v_point_index); + /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if + * there isn't such a point. */ + int pull_near(int u_point_index); + /** \internal \brief Returns and remove all the V points near to the U point given as parameter. */ + std::vector<int> pull_all_near(int u_point_index); + + private: + const Persistence_graph& g; + const double r; + Kd_tree kd_t; + std::unordered_set<int> projections_f; +}; + +/** \internal \brief data structure used to find any point (including projections) in V near to a query point from U + * (which can be a projection) in a layered graph layer given as parmeter. + * + * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically + * removed. + * + * \ingroup bottleneck_distance + */ +class Layered_neighbors_finder { + public: + /** \internal \brief Constructor taking the near distance definition as parameter. */ + Layered_neighbors_finder(const Persistence_graph& g, double r); + /** \internal \brief A point added will be possibly pulled. */ + void add(int v_point_index, int vlayer); + /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if + * there isn't such a point. */ + int pull_near(int u_point_index, int vlayer); + /** \internal \brief Returns the number of layers. */ + int vlayers_number() const; + + private: + const Persistence_graph& g; + const double r; + std::vector<std::unique_ptr<Neighbors_finder>> neighbors_finder; +}; + +inline Neighbors_finder::Neighbors_finder(const Persistence_graph& g, double r) : + g(g), r(r), kd_t(), projections_f() { } + +inline void Neighbors_finder::add(int v_point_index) { + if (g.on_the_v_diagonal(v_point_index)) + projections_f.emplace(v_point_index); + else + kd_t.insert(g.get_v_point(v_point_index)); +} + +inline int Neighbors_finder::pull_near(int u_point_index) { + int tmp; + int c = g.corresponding_point_in_v(u_point_index); + if (g.on_the_u_diagonal(u_point_index) && !projections_f.empty()) { + // Any pair of projection is at distance 0 + tmp = *projections_f.cbegin(); + projections_f.erase(tmp); + } else if (projections_f.count(c) && (g.distance(u_point_index, c) <= r)) { + // Is the query point near to its projection ? + tmp = c; + projections_f.erase(tmp); + } else { + // Is the query point near to a V point in the plane ? + Internal_point u_point = g.get_u_point(u_point_index); + auto neighbor = kd_t.search_any_point(Square_query{u_point, r}); + if (!neighbor) + return null_point_index(); + tmp = neighbor->point_index; + auto point = g.get_v_point(tmp); + int idx = point.point_index; + kd_t.remove(point, [idx](Internal_point const&p){return p.point_index == idx;}); + } + return tmp; +} + +inline std::vector<int> Neighbors_finder::pull_all_near(int u_point_index) { + std::vector<int> all_pull; + int last_pull = pull_near(u_point_index); + while (last_pull != null_point_index()) { + all_pull.emplace_back(last_pull); + last_pull = pull_near(u_point_index); + } + return all_pull; +} + +inline Layered_neighbors_finder::Layered_neighbors_finder(const Persistence_graph& g, double r) : + g(g), r(r), neighbors_finder() { } + +inline void Layered_neighbors_finder::add(int v_point_index, int vlayer) { + for (int l = neighbors_finder.size(); l <= vlayer; l++) + neighbors_finder.emplace_back(std::unique_ptr<Neighbors_finder>(new Neighbors_finder(g, r))); + neighbors_finder.at(vlayer)->add(v_point_index); +} + +inline int Layered_neighbors_finder::pull_near(int u_point_index, int vlayer) { + if (static_cast<int> (neighbors_finder.size()) <= vlayer) + return null_point_index(); + return neighbors_finder.at(vlayer)->pull_near(u_point_index); +} + +inline int Layered_neighbors_finder::vlayers_number() const { + return static_cast<int> (neighbors_finder.size()); +} + +} // namespace persistence_diagram + +} // namespace Gudhi + +#endif // NEIGHBORS_FINDER_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h new file mode 100644 index 00000000..f791e37c --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h @@ -0,0 +1,176 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author: Francois Godi + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef PERSISTENCE_GRAPH_H_ +#define PERSISTENCE_GRAPH_H_ + +#include <gudhi/Internal_point.h> + +#ifdef GUDHI_USE_TBB +#include <tbb/parallel_sort.h> +#endif + +#include <vector> +#include <algorithm> +#include <limits> // for numeric_limits + +namespace Gudhi { + +namespace persistence_diagram { + +/** \internal \brief Structure representing an euclidean bipartite graph containing + * the points from the two persistence diagrams (including the projections). + * + * \ingroup bottleneck_distance + */ +class Persistence_graph { + public: + /** \internal \brief Constructor taking 2 PersistenceDiagrams (concept) as parameters. */ + template<typename Persistence_diagram1, typename Persistence_diagram2> + Persistence_graph(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e); + /** \internal \brief Is the given point from U the projection of a point in V ? */ + bool on_the_u_diagonal(int u_point_index) const; + /** \internal \brief Is the given point from V the projection of a point in U ? */ + bool on_the_v_diagonal(int v_point_index) const; + /** \internal \brief Given a point from V, returns the corresponding (projection or projector) point in U. */ + int corresponding_point_in_u(int v_point_index) const; + /** \internal \brief Given a point from U, returns the corresponding (projection or projector) point in V. */ + int corresponding_point_in_v(int u_point_index) const; + /** \internal \brief Given a point from U and a point from V, returns the distance between those points. */ + double distance(int u_point_index, int v_point_index) const; + /** \internal \brief Returns size = |U| = |V|. */ + int size() const; + /** \internal \brief Is there as many infinite points (alive components) in both diagrams ? */ + double bottleneck_alive() const; + /** \internal \brief Returns the O(n^2) sorted distances between the points. */ + std::vector<double> sorted_distances() const; + /** \internal \brief Returns an upper bound for the diameter of the convex hull of all non infinite points */ + double diameter_bound() const; + /** \internal \brief Returns the corresponding internal point */ + Internal_point get_u_point(int u_point_index) const; + /** \internal \brief Returns the corresponding internal point */ + Internal_point get_v_point(int v_point_index) const; + + private: + std::vector<Internal_point> u; + std::vector<Internal_point> v; + double b_alive; +}; + +template<typename Persistence_diagram1, typename Persistence_diagram2> +Persistence_graph::Persistence_graph(const Persistence_diagram1 &diag1, + const Persistence_diagram2 &diag2, double e) + : u(), v(), b_alive(0.) { + std::vector<double> u_alive; + std::vector<double> v_alive; + for (auto it = std::begin(diag1); it != std::end(diag1); ++it) { + if (std::get<1>(*it) == std::numeric_limits<double>::infinity()) + u_alive.push_back(std::get<0>(*it)); + else if (std::get<1>(*it) - std::get<0>(*it) > e) + u.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), u.size())); + } + for (auto it = std::begin(diag2); it != std::end(diag2); ++it) { + if (std::get<1>(*it) == std::numeric_limits<double>::infinity()) + v_alive.push_back(std::get<0>(*it)); + else if (std::get<1>(*it) - std::get<0>(*it) > e) + v.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), v.size())); + } + if (u.size() < v.size()) + swap(u, v); + std::sort(u_alive.begin(), u_alive.end()); + std::sort(v_alive.begin(), v_alive.end()); + if (u_alive.size() != v_alive.size()) { + b_alive = std::numeric_limits<double>::infinity(); + } else { + for (auto it_u = u_alive.cbegin(), it_v = v_alive.cbegin(); it_u != u_alive.cend(); ++it_u, ++it_v) + b_alive = (std::max)(b_alive, std::fabs(*it_u - *it_v)); + } +} + +inline bool Persistence_graph::on_the_u_diagonal(int u_point_index) const { + return u_point_index >= static_cast<int> (u.size()); +} + +inline bool Persistence_graph::on_the_v_diagonal(int v_point_index) const { + return v_point_index >= static_cast<int> (v.size()); +} + +inline int Persistence_graph::corresponding_point_in_u(int v_point_index) const { + return on_the_v_diagonal(v_point_index) ? + v_point_index - static_cast<int> (v.size()) : v_point_index + static_cast<int> (u.size()); +} + +inline int Persistence_graph::corresponding_point_in_v(int u_point_index) const { + return on_the_u_diagonal(u_point_index) ? + u_point_index - static_cast<int> (u.size()) : u_point_index + static_cast<int> (v.size()); +} + +inline double Persistence_graph::distance(int u_point_index, int v_point_index) const { + if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(v_point_index)) + return 0.; + Internal_point p_u = get_u_point(u_point_index); + Internal_point p_v = get_v_point(v_point_index); + return (std::max)(std::fabs(p_u.x() - p_v.x()), std::fabs(p_u.y() - p_v.y())); +} + +inline int Persistence_graph::size() const { + return static_cast<int> (u.size() + v.size()); +} + +inline double Persistence_graph::bottleneck_alive() const { + return b_alive; +} + +inline std::vector<double> Persistence_graph::sorted_distances() const { + std::vector<double> distances; + distances.push_back(0.); // for empty diagrams + for (int u_point_index = 0; u_point_index < size(); ++u_point_index) { + distances.push_back(distance(u_point_index, corresponding_point_in_v(u_point_index))); + for (int v_point_index = 0; v_point_index < size(); ++v_point_index) + distances.push_back(distance(u_point_index, v_point_index)); + } +#ifdef GUDHI_USE_TBB + tbb::parallel_sort(distances.begin(), distances.end()); +#else + std::sort(distances.begin(), distances.end()); +#endif + return distances; +} + +inline Internal_point Persistence_graph::get_u_point(int u_point_index) const { + if (!on_the_u_diagonal(u_point_index)) + return u.at(u_point_index); + Internal_point projector = v.at(corresponding_point_in_v(u_point_index)); + double m = (projector.x() + projector.y()) / 2.; + return Internal_point(m, m, u_point_index); +} + +inline Internal_point Persistence_graph::get_v_point(int v_point_index) const { + if (!on_the_v_diagonal(v_point_index)) + return v.at(v_point_index); + Internal_point projector = u.at(corresponding_point_in_u(v_point_index)); + double m = (projector.x() + projector.y()) / 2.; + return Internal_point(m, m, v_point_index); +} + +inline double Persistence_graph::diameter_bound() const { + double max = 0.; + for (auto it = u.cbegin(); it != u.cend(); it++) + max = (std::max)(max, it->y()); + for (auto it = v.cbegin(); it != v.cend(); it++) + max = (std::max)(max, it->y()); + return max; +} + +} // namespace persistence_diagram + +} // namespace Gudhi + +#endif // PERSISTENCE_GRAPH_H_ diff --git a/src/Bottleneck_distance/test/CMakeLists.txt b/src/Bottleneck_distance/test/CMakeLists.txt new file mode 100644 index 00000000..ec2d045f --- /dev/null +++ b/src/Bottleneck_distance/test/CMakeLists.txt @@ -0,0 +1,14 @@ +project(Bottleneck_distance_tests) + +if (NOT CGAL_VERSION VERSION_LESS 4.11.0) + include(GUDHI_test_coverage) + + add_executable ( Bottleneck_distance_test_unit bottleneck_unit_test.cpp ) + target_link_libraries(Bottleneck_distance_test_unit ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + if (TBB_FOUND) + target_link_libraries(Bottleneck_distance_test_unit ${TBB_LIBRARIES}) + endif(TBB_FOUND) + + gudhi_add_coverage_test(Bottleneck_distance_test_unit) + +endif (NOT CGAL_VERSION VERSION_LESS 4.11.0) diff --git a/src/Bottleneck_distance/test/README b/src/Bottleneck_distance/test/README new file mode 100644 index 00000000..00bd9a91 --- /dev/null +++ b/src/Bottleneck_distance/test/README @@ -0,0 +1,12 @@ +To compile: +*********** + +cmake . +make + +To launch with details: +*********************** + +./Bottleneck_distance_unit_test --report_level=detailed --log_level=all + + ==> echo $? returns 0 in case of success (non-zero otherwise) diff --git a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp new file mode 100644 index 00000000..3fc6fc7b --- /dev/null +++ b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp @@ -0,0 +1,155 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author: Francois Godi + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "bottleneck distance" +#include <boost/test/unit_test.hpp> + +#include <random> +#include <gudhi/Bottleneck.h> + +using namespace Gudhi::persistence_diagram; + +int n1 = 81; // a natural number >0 +int n2 = 180; // a natural number >0 +double upper_bound = 406.43; // any real >0 + + +std::uniform_real_distribution<double> unif(0., upper_bound); +std::default_random_engine re; +std::vector< std::pair<double, double> > v1, v2; + +BOOST_AUTO_TEST_CASE(persistence_graph) { + // Random construction + for (int i = 0; i < n1; i++) { + double a = unif(re); + double b = unif(re); + v1.emplace_back(std::min(a, b), std::max(a, b)); + } + for (int i = 0; i < n2; i++) { + double a = unif(re); + double b = unif(re); + v2.emplace_back(std::min(a, b), std::max(a, b)); + } + Persistence_graph g(v1, v2, 0.); + std::vector<double> d(g.sorted_distances()); + // + BOOST_CHECK(!g.on_the_u_diagonal(n1 - 1)); + BOOST_CHECK(!g.on_the_u_diagonal(n1)); + BOOST_CHECK(!g.on_the_u_diagonal(n2 - 1)); + BOOST_CHECK(g.on_the_u_diagonal(n2)); + BOOST_CHECK(!g.on_the_v_diagonal(n1 - 1)); + BOOST_CHECK(g.on_the_v_diagonal(n1)); + BOOST_CHECK(g.on_the_v_diagonal(n2 - 1)); + BOOST_CHECK(g.on_the_v_diagonal(n2)); + // + BOOST_CHECK(g.corresponding_point_in_u(0) == n2); + BOOST_CHECK(g.corresponding_point_in_u(n1) == 0); + BOOST_CHECK(g.corresponding_point_in_v(0) == n1); + BOOST_CHECK(g.corresponding_point_in_v(n2) == 0); + // + BOOST_CHECK(g.size() == (n1 + n2)); + // + BOOST_CHECK((int) d.size() == (n1 + n2)*(n1 + n2) + n1 + n2 + 1); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, 0)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, n1 - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, n1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, n2 - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, n2)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, (n1 + n2) - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, 0)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, n1 - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, n1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, n2 - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, n2)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, (n1 + n2) - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, 0)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, n1 - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, n1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, n2 - 1)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, n2)) > 0); + BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, (n1 + n2) - 1)) > 0); +} + +BOOST_AUTO_TEST_CASE(neighbors_finder) { + Persistence_graph g(v1, v2, 0.); + Neighbors_finder nf(g, 1.); + for (int v_point_index = 1; v_point_index < ((n2 + n1)*9 / 10); v_point_index += 2) + nf.add(v_point_index); + // + int v_point_index_1 = nf.pull_near(n2 / 2); + BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2 / 2, v_point_index_1) <= 1.)); + std::vector<int> l = nf.pull_all_near(n2 / 2); + bool v = true; + for (auto it = l.cbegin(); it != l.cend(); ++it) + v = v && (g.distance(n2 / 2, *it) > 1.); + BOOST_CHECK(v); + int v_point_index_2 = nf.pull_near(n2 / 2); + BOOST_CHECK(v_point_index_2 == -1); +} + +BOOST_AUTO_TEST_CASE(layered_neighbors_finder) { + Persistence_graph g(v1, v2, 0.); + Layered_neighbors_finder lnf(g, 1.); + for (int v_point_index = 1; v_point_index < ((n2 + n1)*9 / 10); v_point_index += 2) + lnf.add(v_point_index, v_point_index % 7); + // + int v_point_index_1 = lnf.pull_near(n2 / 2, 6); + BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2 / 2, v_point_index_1) <= 1.)); + int v_point_index_2 = lnf.pull_near(n2 / 2, 6); + BOOST_CHECK(v_point_index_2 == -1); + v_point_index_1 = lnf.pull_near(n2 / 2, 0); + BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2 / 2, v_point_index_1) <= 1.)); + v_point_index_2 = lnf.pull_near(n2 / 2, 0); + BOOST_CHECK(v_point_index_2 == -1); +} + +BOOST_AUTO_TEST_CASE(graph_matching) { + Persistence_graph g(v1, v2, 0.); + Graph_matching m1(g); + m1.set_r(0.); + int e = 0; + while (m1.multi_augment()) + ++e; + BOOST_CHECK(e > 0); + BOOST_CHECK(e <= 2 * sqrt(2 * (n1 + n2))); + Graph_matching m2 = m1; + BOOST_CHECK(!m2.multi_augment()); + m2.set_r(upper_bound); + e = 0; + while (m2.multi_augment()) + ++e; + BOOST_CHECK(e <= 2 * sqrt(2 * (n1 + n2))); + BOOST_CHECK(m2.perfect()); + BOOST_CHECK(!m1.perfect()); +} + +BOOST_AUTO_TEST_CASE(global) { + std::uniform_real_distribution<double> unif1(0., upper_bound); + std::uniform_real_distribution<double> unif2(upper_bound / 10000., upper_bound / 100.); + std::default_random_engine re; + std::vector< std::pair<double, double> > v1, v2; + for (int i = 0; i < n1; i++) { + double a = unif1(re); + double b = unif1(re); + double x = unif2(re); + double y = unif2(re); + v1.emplace_back(std::min(a, b), std::max(a, b)); + v2.emplace_back(std::min(a, b) + std::min(x, y), std::max(a, b) + std::max(x, y)); + if (i % 5 == 0) + v1.emplace_back(std::min(a, b), std::min(a, b) + x); + if (i % 3 == 0) + v2.emplace_back(std::max(a, b), std::max(a, b) + y); + } + BOOST_CHECK(bottleneck_distance(v1, v2, 0.) <= upper_bound / 100.); + BOOST_CHECK(bottleneck_distance(v1, v2, upper_bound / 10000.) <= upper_bound / 100. + upper_bound / 10000.); + BOOST_CHECK(std::abs(bottleneck_distance(v1, v2, 0.) - bottleneck_distance(v1, v2, upper_bound / 10000.)) <= upper_bound / 10000.); +} diff --git a/src/Bottleneck_distance/utilities/CMakeLists.txt b/src/Bottleneck_distance/utilities/CMakeLists.txt new file mode 100644 index 00000000..86d74cf5 --- /dev/null +++ b/src/Bottleneck_distance/utilities/CMakeLists.txt @@ -0,0 +1,15 @@ +project(Bottleneck_distance_utilities) + +if(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + add_executable (bottleneck_distance bottleneck_distance.cpp) + if (TBB_FOUND) + target_link_libraries(bottleneck_distance ${TBB_LIBRARIES}) + endif(TBB_FOUND) + + add_test(NAME Bottleneck_distance_utilities_Bottleneck_read_file + COMMAND $<TARGET_FILE:bottleneck_distance> + "${CMAKE_SOURCE_DIR}/data/persistence_diagram/first.pers" "${CMAKE_SOURCE_DIR}/data/persistence_diagram/second.pers") + + install(TARGETS bottleneck_distance DESTINATION bin) + +endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) diff --git a/src/Bottleneck_distance/utilities/bottleneck_distance.cpp b/src/Bottleneck_distance/utilities/bottleneck_distance.cpp new file mode 100644 index 00000000..d88a8a0b --- /dev/null +++ b/src/Bottleneck_distance/utilities/bottleneck_distance.cpp @@ -0,0 +1,38 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Authors: Francois Godi, small modifications by Pawel Dlotko + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <gudhi/Bottleneck.h> +#include <gudhi/reader_utils.h> +#include <iostream> +#include <vector> +#include <utility> // for pair +#include <string> +#include <limits> // for numeric_limits + +int main(int argc, char** argv) { + if (argc < 3) { + std::cout << "To run this program please provide as an input two files with persistence diagrams. Each file" << + " should contain a birth-death pair per line. Third, optional parameter is an error bound on the bottleneck" << + " distance (set by default to the smallest positive double value). If you set the error bound to 0, be" << + " aware this version is exact but expensive. The program will now terminate \n"; + return -1; + } + std::vector<std::pair<double, double>> diag1 = Gudhi::read_persistence_intervals_in_dimension(argv[1]); + std::vector<std::pair<double, double>> diag2 = Gudhi::read_persistence_intervals_in_dimension(argv[2]); + + double tolerance = std::numeric_limits<double>::min(); + if (argc == 4) { + tolerance = atof(argv[3]); + } + double b = Gudhi::persistence_diagram::bottleneck_distance(diag1, diag2, tolerance); + std::cout << "The distance between the diagrams is : " << b << ". The tolerance is : " << tolerance << std::endl; + + return 0; +} diff --git a/src/Bottleneck_distance/utilities/bottleneckdistance.md b/src/Bottleneck_distance/utilities/bottleneckdistance.md new file mode 100644 index 00000000..a81426cf --- /dev/null +++ b/src/Bottleneck_distance/utilities/bottleneckdistance.md @@ -0,0 +1,26 @@ +---
+layout: page
+title: "Bottleneck distance"
+meta_title: "Bottleneck distance"
+teaser: ""
+permalink: /bottleneckdistance/
+---
+{::comment}
+Leave the lines above as it is required by the web site generator 'Jekyll'
+{:/comment}
+
+
+## bottleneck_read_file_example ##
+
+This program computes the Bottleneck distance between two persistence diagram files.
+
+**Usage**
+
+```
+ bottleneck_read_file_example <file_1.pers> <file_2.pers> [<tolerance>]
+```
+
+where
+
+* `<file_1.pers>` and `<file_2.pers>` must be in the format described [here]({{ site.officialurl }}/doc/latest/fileformats.html#FileFormatsPers).
+* `<tolerance>` is an error bound on the bottleneck distance (set by default to the smallest positive double value).
|
