diff --git a/CMakeLists.txt b/CMakeLists.txt index 02cdbe61c9..a0752b637a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -47,6 +47,7 @@ add_gudhi_module(Toplex_map) add_gudhi_module(Witness_complex) add_gudhi_module(Nerve_GIC) add_gudhi_module(Persistence_matrix) +add_gudhi_module(Zigzag_persistence) # Include module CMake subdirectories # GUDHI_SUB_DIRECTORIES is managed in CMAKE_MODULE_PATH/GUDHI_modules.cmake diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index bae5a0ee7f..bed9459fc2 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -46,6 +46,7 @@ add_gudhi_module(Toplex_map) add_gudhi_module(Witness_complex) add_gudhi_module(Nerve_GIC) add_gudhi_module(Persistence_matrix) +add_gudhi_module(Zigzag_persistence) set(GUDHI_BIBLIO_DIR ${CMAKE_SOURCE_DIR}) # For "make doxygen" - Requires GUDHI_USER_VERSION_DIR to be set diff --git a/src/Zigzag_persistence/concept/ZigzagOptions.h b/src/Zigzag_persistence/concept/ZigzagOptions.h new file mode 100644 index 0000000000..35f421ac11 --- /dev/null +++ b/src/Zigzag_persistence/concept/ZigzagOptions.h @@ -0,0 +1,80 @@ +/* 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): Hannah Schreiber + * + * Copyright (C) 2024 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef CONCEPT_ZZ_OPTIONS_TYPE_H_ +#define CONCEPT_ZZ_OPTIONS_TYPE_H_ + +/** @file ZigzagOptions.h + * @brief Contains @ref Gudhi::zigzag_persistence::ZigzagOptions and + * @ref Gudhi::zigzag_persistence::FilteredZigzagOptions concept. + */ + +namespace Gudhi { +namespace zigzag_persistence { + +/** + * @ingroup zigzag_persistence + * + * @brief List of options used for the filtered zigzag persistence computation. + */ +struct FilteredZigzagOptions { + /** + * @brief Numerical type for the cell IDs used internally and other indexations. It must be signed. + */ + using Internal_key = unspecified; + + /** + * @brief Type for the cell IDs used at insertion and in the boundaries given as argument. + * Has to be usable as key in a hashtable, so "hashable" and comparable. + */ + using Cell_key = unspecified; + + /** + * @brief Type for filtration values. + */ + using Filtration_value = unspecified; + + /** + * @brief Type for the dimension values. + */ + using Dimension = unspecified; + + /** + * @brief Column type used by the internal matrix. + */ + static const Gudhi::persistence_matrix::Column_types column_type; +}; + +/** + * @ingroup zigzag_persistence + * + * @brief List of options used for the zigzag persistence computation. + */ +struct ZigzagOptions { + /** + * @brief Numerical type for the cell IDs used internally and other indexations. It must be signed. + */ + using Internal_key = unspecified; + + /** + * @brief Type for the dimension values. + */ + using Dimension = unspecified; + + /** + * @brief Column type used by the internal matrix. + */ + static const Gudhi::persistence_matrix::Column_types column_type; +}; + +} // namespace zigzag_persistence +} // namespace Gudhi + +#endif // CONCEPT_ZZ_OPTIONS_TYPE_H_ diff --git a/src/Zigzag_persistence/doc/COPYRIGHT b/src/Zigzag_persistence/doc/COPYRIGHT new file mode 100644 index 0000000000..f9ac497dac --- /dev/null +++ b/src/Zigzag_persistence/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): Hannah Schreiber + +Copyright (C) 2024 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/Zigzag_persistence/doc/Intro_zigzag_persistence.h b/src/Zigzag_persistence/doc/Intro_zigzag_persistence.h new file mode 100644 index 0000000000..282825e7ba --- /dev/null +++ b/src/Zigzag_persistence/doc/Intro_zigzag_persistence.h @@ -0,0 +1,97 @@ +/* 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): Hannah Schreiber + * + * Copyright (C) 2023 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef DOC_ZIGZAG_PERSISTENCE_INTRO_ZIGZAG_PERSISTENCE_H_ +#define DOC_ZIGZAG_PERSISTENCE_INTRO_ZIGZAG_PERSISTENCE_H_ + +// needs namespace for Doxygen to link on classes +namespace Gudhi { +namespace zigzag_persistence { + +/** \defgroup zigzag_persistence Zigzag Persistence + * @{ + * \author Clément Maria, Hannah Schreiber + * + * \section zigzagintro Zigzag Persistence + * + * We refer to the introduction page \ref persistent_cohomology for persistent (co)homology for an introduction + * to the topic. + * Zigzag persistence is a generalization of the latter. While standard persistence only allows to grow the filtered + * complex by adding cells, zigzag persistence also allows removals. Hence the name "zigzag", as the module + * diagram will have arrows alternating between forward and backward. + * + * The module consists of the @ref Zigzag_persistence class and two wrappers @ref Filtered_zigzag_persistence and + * @ref Filtered_zigzag_persistence_with_storage "": + * - @ref Zigzag_persistence computes the persistence of a sequence of insertions and removals. A cell can be inserted + * or removed one at a time and the returned persistence pairs / bars are indexed on the operation numbers. + * For example, if a cycle is born at operation number 6 and dies at operation number 7, it will output a bar starting + * at 6 and ending at 7. + * - @ref Filtered_zigzag_persistence and @ref Filtered_zigzag_persistence_with_storage are adding the notion of + * "filtration value" to @ref Zigzag_persistence. At each call, an operation can be associated to a filtration value, + * which will be used to index the returned bars instead (bars with new length 0 are then ignored). The two classes + * also have more flexible inputs (the boundaries do not have to be ordered, nor identified continuously + * from 0). The difference between both classes is on the way they manage the memory: @ref Filtered_zigzag_persistence + * removes systematically all unnecessary information and outputs a pair as soon it is closed, while + * @ref Filtered_zigzag_persistence_with_storage will store all informations about filtration values and bars until the + * end and output the pairs only when asked. Depending on the use and the length of the filtration, one will be more + * efficient than the other and vice versa. + * + * The implementation is based on the algorithm introduced in \cite zigzag. + * + * \subsection zigzaginterface Stream-like interface + * + * As removals are possible in zigzag filtration, the maximal size of the complex does not depend on the length of the + * filtration anymore. This makes it possible to build very long fine tuned filtrations with relatively small complexes + * which can be processed without overreaching memory space. For this purpose, it is possible to feed the module with + * information about the filtration "on the fly" to avoid loading the whole filtration at once. Information about the + * current barcode can be retrieved between any steps via callback methods. + * + * \section zigzagexamples Examples + * + * \subsection zzminusage Minimalistic examples + * + * \li \gudhi_example_link{Zigzag_persistence,example_usage_zigzag_persistence.cpp} - A simple example to showcase how + * to use the @ref Zigzag_persistence class to compute a barcode. + *
+ * @dontinclude example_usage_zigzag_persistence.cpp + * @skip #include + * @until return 0; + * @skipline } + *
+ * \li \gudhi_example_link{Zigzag_persistence,example_usage_filtered_zigzag_persistence.cpp} - A simple example to + * showcase how to use the @ref Filtered_zigzag_persistence class to compute a barcode. + *
+ * @dontinclude example_usage_filtered_zigzag_persistence.cpp + * @skip #include + * @until return 0; + * @skipline } + *
+ * \li \gudhi_example_link{Zigzag_persistence,example_usage_filtered_zigzag_persistence_with_storage.cpp} - A simple + * example to showcase how to use the @ref Filtered_zigzag_persistence_with_storage class to compute a barcode. + *
+ * @dontinclude example_usage_filtered_zigzag_persistence_with_storage.cpp + * @skip #include + * @until return 0; + * @skipline } + *
+ * + * \subsection zzexamples More elaborate examples + * + * \li \gudhi_example_link{Zigzag_persistence,example_zigzag_filtration_as_input_loop.cpp} - A simple example to showcase how + * to use the @ref Filtered_zigzag_persistence_with_storage class within an input loop. + * \li \gudhi_example_link{Zigzag_persistence,example_zzfiltration_from_file.cpp} - An example of a "stream-like" usage + * with @ref Filtered_zigzag_persistence by reading off the filtration from a file. + * + * @} + */ +} // namespace zigzag_persistence +} // namespace Gudhi + +#endif // DOC_ZIGZAG_PERSISTENCE_INTRO_ZIGZAG_PERSISTENCE_H_ diff --git a/src/Zigzag_persistence/doc/zigzag_ex.png b/src/Zigzag_persistence/doc/zigzag_ex.png new file mode 100644 index 0000000000..ba7cb7011a Binary files /dev/null and b/src/Zigzag_persistence/doc/zigzag_ex.png differ diff --git a/src/Zigzag_persistence/example/CMakeLists.txt b/src/Zigzag_persistence/example/CMakeLists.txt new file mode 100644 index 0000000000..ac434df9e6 --- /dev/null +++ b/src/Zigzag_persistence/example/CMakeLists.txt @@ -0,0 +1,21 @@ +project(Zigzag_persistence_examples) + +add_executable_with_targets(Zigzag_persistence_example_usage_zigzag_persistence example_usage_zigzag_persistence.cpp TBB::tbb) +add_test(NAME Zigzag_persistence_example_usage_zigzag_persistence COMMAND $) + +add_executable_with_targets(Zigzag_persistence_example_usage_filtered_zigzag_persistence example_usage_filtered_zigzag_persistence.cpp TBB::tbb) +add_test(NAME Zigzag_persistence_example_usage_filtered_zigzag_persistence COMMAND $) + +add_executable_with_targets(Zigzag_persistence_example_usage_filtered_zigzag_persistence_with_storage example_usage_filtered_zigzag_persistence_with_storage.cpp TBB::tbb) +add_test(NAME Zigzag_persistence_example_usage_filtered_zigzag_persistence_with_storage COMMAND $) + +add_executable_with_targets(Zigzag_persistence_example_zigzag_filtration_as_input_loop example_zigzag_filtration_as_input_loop.cpp TBB::tbb) +add_test(NAME Zigzag_persistence_example_zigzag_filtration_as_input_loop COMMAND $) + +add_executable_with_targets(Zigzag_persistence_example_zzfiltration_from_file example_zzfiltration_from_file.cpp TBB::tbb) +file(COPY "zigzag_filtration_example.txt" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) +add_test(NAME Zigzag_persistence_example_zzfiltration_from_file COMMAND $ "${CMAKE_CURRENT_BINARY_DIR}/zigzag_filtration_example.txt") + + + + diff --git a/src/Zigzag_persistence/example/example_usage_filtered_zigzag_persistence.cpp b/src/Zigzag_persistence/example/example_usage_filtered_zigzag_persistence.cpp new file mode 100644 index 0000000000..1ceb362015 --- /dev/null +++ b/src/Zigzag_persistence/example/example_usage_filtered_zigzag_persistence.cpp @@ -0,0 +1,57 @@ +/* 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): Hannah Schreiber + * + * Copyright (C) 2024 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include + +#include + +using Zigzag_persistence = Gudhi::zigzag_persistence::Filtered_zigzag_persistence<>; +using Dimension = Zigzag_persistence::Dimension; +using Filtration_value = Zigzag_persistence::Filtration_value; + +int main() { + std::clog << "********* Minimalistic example of usage of the Filtered_zigzag_persistence class ********" << std::endl; + + // Filtered_zigzag_persistence(callback) with for example callback method as a anonymous lambda + Zigzag_persistence zp([](Dimension dim, Filtration_value birth, Filtration_value death) { + std::cout << "[" << dim << "] " << birth << " - " << death << std::endl; + }); + + // It is important that the operations of insertions and removals are made **in the same order** as in the zigzag + // filtration ones wants to compute the barcode from. + // A cell can be identified in the boundaries by any given numerical label, it is just important that the given + // filtration values are monotonous (ie., either only increasing or only decreasing). + + // inserts vertex 2 at filtration value 0.1 -> birth at 0.1 of 0-cycle + zp.insert_cell(2, {}, 0, 0.1); + // inserts vertex 4 at filtration value 0.1 -> birth at 0.1 of 0-cycle + zp.insert_cell(4, {}, 0, 0.1); + // inserts edge 5 = (2,4) at filtration value 0.3 -> death at 0.3 -> outputs (0, 0.1, 0.3) + zp.insert_cell(5, {2, 4}, 1, 0.3); + // inserts vertex 3 at filtration value 0.4 -> birth at 0.4 of 0-cycle + zp.insert_cell(3, {}, 0, 0.4); + // inserts edge 6 = (2,3) at filtration value 0.4 -> death at 0.4 of the cycle born at 0.4 -> outputs nothing + zp.insert_cell(6, {2, 3}, 1, 0.4); + // inserts edge 9 = (3,4) at filtration value 1.2 -> birth at 1.2 of 1-cycle + zp.insert_cell(9, {4, 3}, 1, 1.2); + // removes edge 6 at filtration value 1.5 -> death at 1.5 -> outputs (1, 1.2, 1.5) + zp.remove_cell(6, 1.5); + // removes edge 5 at filtration value 2.0 -> birth at 2.0 of 0-cycle + zp.remove_cell(5, 2.0); + + // Only the closed bars where output so far, so the open/infinite bars still need to be retrieved. + + // in this example, computes (0, 0.1) and (0, 2.0) + zp.get_current_infinite_intervals([](Dimension dim, Filtration_value birth) { + std::cout << "[" << dim << "] " << birth << " - inf" << std::endl; + }); + + return 0; +} diff --git a/src/Zigzag_persistence/example/example_usage_filtered_zigzag_persistence_with_storage.cpp b/src/Zigzag_persistence/example/example_usage_filtered_zigzag_persistence_with_storage.cpp new file mode 100644 index 0000000000..cd14b60fe0 --- /dev/null +++ b/src/Zigzag_persistence/example/example_usage_filtered_zigzag_persistence_with_storage.cpp @@ -0,0 +1,55 @@ +/* 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): Hannah Schreiber + * + * Copyright (C) 2024 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include + +#include + +using Zigzag_persistence = Gudhi::zigzag_persistence::Filtered_zigzag_persistence_with_storage<>; + +int main() { + std::clog << "** Minimalistic example of usage of the Filtered_zigzag_persistence_with_storage class **" << std::endl; + + Zigzag_persistence zp; + + // It is important that the operations of insertions and removals are made **in the same order** as in the zigzag + // filtration ones wants to compute the barcode from. + // A cell can be identified in the boundaries by any given numerical label, it is just important that the given + // filtration values are monotonous (ie., either only increasing or only decreasing). + + // inserts vertex 2 at filtration value 0.1 -> birth at 0.1 of 0-cycle + zp.insert_cell(2, {}, 0, 0.1); + // inserts vertex 4 at filtration value 0.1 -> birth at 0.1 of 0-cycle + zp.insert_cell(4, {}, 0, 0.1); + // inserts edge 5 = (2,4) at filtration value 0.3 -> death at 0.3 -> stores (0, 0.1, 0.3) + zp.insert_cell(5, {2, 4}, 1, 0.3); + // inserts vertex 3 at filtration value 0.4 -> birth at 0.4 of 0-cycle + zp.insert_cell(3, {}, 0, 0.4); + // inserts edge 6 = (2,3) at filtration value 0.4 -> death at 0.4 of the cycle born at 0.4 -> stores nothing + zp.insert_cell(6, {2, 3}, 1, 0.4); + // inserts edge 9 = (3,4) at filtration value 1.2 -> birth at 1.2 of 1-cycle + zp.insert_cell(9, {4, 3}, 1, 1.2); + // removes edge 6 at filtration value 1.5 -> death at 1.5 -> stores (1, 1.2, 1.5) + zp.remove_cell(6, 1.5); + // removes edge 5 at filtration value 2.0 -> birth at 2.0 of 0-cycle + zp.remove_cell(5, 2.0); + + // The bars are stored within the class and where not output at all for now. + + // get all bars in a vector + auto barcode = zp.get_persistence_diagram(); + + // do something with the vector, e.g., stream out content: + for (auto& bar : barcode) { + std::cout << bar << std::endl; + } + + return 0; +} diff --git a/src/Zigzag_persistence/example/example_usage_zigzag_persistence.cpp b/src/Zigzag_persistence/example/example_usage_zigzag_persistence.cpp new file mode 100644 index 0000000000..e50f94cfc0 --- /dev/null +++ b/src/Zigzag_persistence/example/example_usage_zigzag_persistence.cpp @@ -0,0 +1,55 @@ +/* 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): Hannah Schreiber + * + * Copyright (C) 2024 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include + +#include + +using Zigzag_persistence = Gudhi::zigzag_persistence::Zigzag_persistence<>; +using Dimension = Zigzag_persistence::Dimension; +using Index = Zigzag_persistence::Index; + +int main() { + std::clog << "************* Minimalistic example of usage of the Zigzag_persistence class *************" << std::endl; + + // Zigzag_persistence(callback) with for example callback method as a anonymous lambda + Zigzag_persistence zp([](Dimension dim, Index birth, Index death) { + std::cout << "[" << dim << "] " << birth << " - " << death << std::endl; + }); + + // It is important that the operations of insertions and removals are made **in the same order** as in the zigzag + // filtration ones wants to compute the barcode from. + // A cell has to be identified in the boundaries by the operation number the cell was inserted with in the sequence. + + // inserts vertex 0 -> birth at 0 of 0-cycle + zp.insert_cell({}, 0); + // inserts vertex 1 -> birth at 1 of 0-cycle + zp.insert_cell({}, 0); + // inserts edge 2 = (0,1) -> death at 2 -> outputs (0, 1, 2) + zp.insert_cell({0, 1}, 1); + // inserts vertex 3 -> birth at 3 of 0-cycle + zp.insert_cell({}, 0); + // inserts edge 4 = (0,3) -> death at 4 -> outputs (0, 3, 4) + zp.insert_cell({0, 3}, 1); + // inserts edge 5 = (1,3) -> birth at 5 of 1-cycle + zp.insert_cell({1, 3}, 1); + // removes edge 4 -> death at 6 -> outputs (1, 5, 6) + zp.remove_cell(4); + // removes edge 2 -> birth at 7 of 0-cycle + zp.remove_cell(2); + + // Only the closed bars were output so far, so the open/infinite bars still need to be retrieved. + + // in this example, computes (0, 0) and (0, 7) + zp.get_current_infinite_intervals( + [](Dimension dim, Index birth) { std::cout << "[" << dim << "] " << birth << " - inf" << std::endl; }); + + return 0; +} diff --git a/src/Zigzag_persistence/example/example_zigzag_filtration_as_input_loop.cpp b/src/Zigzag_persistence/example/example_zigzag_filtration_as_input_loop.cpp new file mode 100644 index 0000000000..e98dae969e --- /dev/null +++ b/src/Zigzag_persistence/example/example_zigzag_filtration_as_input_loop.cpp @@ -0,0 +1,139 @@ +/* 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): Hannah Schreiber + * + * Copyright (C) 2023 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include +#include + +#include + +using ZP = Gudhi::zigzag_persistence::Filtered_zigzag_persistence_with_storage<>; +using Cell_handle = ZP::Cell_key; +using Filtration_value = ZP::Filtration_value; +using Interval_filtration = ZP::Filtration_value_interval; + +void print_barcode(ZP& zp) { + std::cout << std::endl << "Current barcode:" << std::endl; + for (Interval_filtration& bar : zp.get_persistence_diagram(0, true)) { + //stream out content of bar + std::cout << bar << std::endl; + //to access the content of the bar, it can either be used as a struct: + // bar.birth + // bar.death + // bar.dim + //or as a tuple + // std::get<0>(bar) <- birth + // std::get<1>(bar) <- death + // std::get<2>(bar) <- dim + } +} + +void print_indices(ZP& zp) { + std::cout << std::endl << "Current pairs:" << std::endl; + for (auto& bar : zp.get_index_persistence_diagram()) { + //stream out content of bar + std::cout << bar << std::endl; + //to access the content of the bar, it can either be used as a struct: + // bar.birth + // bar.death + // bar.dim + //or as a tuple: + // std::get<0>(bar) <- birth + // std::get<1>(bar) <- death + // std::get<2>(bar) <- dim + } +} + +std::vector > get_boundaries() { + return {{}, + {}, + {}, + {0, 1}, + {0, 2}, + {}, + {1, 2}, + {}, + {5, 7}, + {}, + {3, 4, 6}, + {7, 9}, + {5, 9}, + {8, 11, 12}, + {10}, // remove + {13}, // remove + {1, 7}, + {3, 4, 6}, + {2, 7}, + {8, 11, 12}, + {0, 7}, + {4, 18, 20}, + {6, 16, 18}, + {3, 16, 20}, + {19}, // remove + {8}, // remove + {12}, // remove + {17, 21, 22, 23}, + {27}}; // remove +} + +std::vector get_filtration_values() { + return {0, 0, 0, + 1, 1, 1, + 2, 2, 2, + 3, 3, 3, 3, + 4, + 5, + 6, 6, 6, + 7, 7, 7, 7, 7, 7, + 8, + 9, 9, 9, + 10}; +} + +std::vector get_directions() { + return {true, true, true, true, true, true, true, true, true, true, true, true, true, true, + false, false, + true, true, true, true, true, true, true, true, + false, false, false, + true, + false}; +} + +std::vector get_batch_sizes() { + return {14, 2, 8, 3, 1, 1}; +} + +int main(int argc, char* const argv[]) { + std::clog << "********** Example **********" << std::endl; + + ZP zp; + + std::vector > simplices = get_boundaries(); + std::vector fils = get_filtration_values(); + std::vector dirs = get_directions(); + + for (unsigned int i = 0; i < simplices.size(); ++i) { + if (i > 0 && dirs[i] != dirs[i - 1]) { + print_barcode(zp); + print_indices(zp); + } + if (dirs[i]) { + int dim = simplices[i].size() == 0 ? 0 : simplices[i].size() - 1; + zp.insert_cell(i, simplices[i], dim, fils[i]); + } else { + auto id = simplices[i][0]; + zp.remove_cell(id, fils[i]); + } + } + + print_barcode(zp); + print_indices(zp); + + return 0; +} diff --git a/src/Zigzag_persistence/example/example_zzfiltration_from_file.cpp b/src/Zigzag_persistence/example/example_zzfiltration_from_file.cpp new file mode 100644 index 0000000000..a3aff58d68 --- /dev/null +++ b/src/Zigzag_persistence/example/example_zzfiltration_from_file.cpp @@ -0,0 +1,125 @@ +/* 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): Hannah Schreiber + * + * Copyright (C) 2023 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include +#include +#include +#include + +#include + +using ZP = Gudhi::zigzag_persistence::Filtered_zigzag_persistence<>; +using ID_handle = ZP::Cell_key; +using Filtration_value = ZP::Filtration_value; +using Dimension = ZP::Dimension; + +enum lineType : int { INCLUSION, REMOVAL, COMMENT }; + +lineType read_operation(std::string& line, std::vector& cells, double& timestamp) { + lineType type; + cells.clear(); + ID_handle num; + + size_t current = line.find_first_not_of(' ', 0); + if (current == std::string::npos) return COMMENT; + + if (line[current] == 'i') + type = INCLUSION; + else if (line[current] == 'r') + type = REMOVAL; + else if (line[current] == '#') + return COMMENT; + else { + std::clog << "(1) Syntaxe error in file." << std::endl; + exit(0); + } + + current = line.find_first_not_of(' ', current + 1); + if (current == std::string::npos) { + std::clog << "(2) Syntaxe error in file." << std::endl; + exit(0); + } + size_t next = line.find_first_of(' ', current); + timestamp = std::stod(line.substr(current, next - current)); + + current = line.find_first_not_of(' ', next); + while (current != std::string::npos) { + next = line.find_first_of(' ', current); + num = std::stoi(line.substr(current, next - current)); + cells.push_back(num); + current = line.find_first_not_of(' ', next); + } + + return type; +} + +//example of input file: example/zigzag_filtration_example.txt +int main(int argc, char* const argv[]) { + if (argc != 2) { + if (argc < 2) + std::clog << "Missing argument: input file name is needed." << std::endl; + else + std::clog << "Too many arguments: only input file name is needed." << std::endl; + return 0; + } + + std::string line; + std::ifstream file(argv[1]); + + //std::cout could be replaced by any other output stream + ZP zp([](Dimension dim, Filtration_value birth, Filtration_value death) { + std::cout << "[" << dim << "] "; + std::cout << birth << " - " << death; + std::cout << std::endl; + }); + + if (file.is_open()) { + std::vector data; + unsigned int id = 0; + double timestamp; + lineType type; + + while (getline(file, line, '\n') && read_operation(line, data, timestamp) == COMMENT); + double lastTimestamp = timestamp; + // first operation has to be an insertion. + zp.insert_cell(id, data, 0, timestamp); + + while (getline(file, line, '\n')) { + type = read_operation(line, data, timestamp); + if (type != COMMENT && lastTimestamp != timestamp) { + lastTimestamp = timestamp; + } + + if (type == INCLUSION) { + ++id; + int dim = data.size() == 0 ? 0 : data.size() - 1; + zp.insert_cell(id, data, dim, timestamp); + } else if (type == REMOVAL) { + ++id; + zp.remove_cell(data[0], timestamp); + } + } + + file.close(); + } else { + std::clog << "Unable to open input file." << std::endl; + file.setstate(std::ios::failbit); + } + + //retrieve infinite bars remaining at the end + //again std::cout could be replaced by any other output stream + zp.get_current_infinite_intervals([](Dimension dim, Filtration_value birth) { + std::cout << "[" << dim << "] "; + std::cout << birth << " - inf"; + std::cout << std::endl; + }); + + return 0; +} \ No newline at end of file diff --git a/src/Zigzag_persistence/example/zigzag_filtration_example.txt b/src/Zigzag_persistence/example/zigzag_filtration_example.txt new file mode 100644 index 0000000000..e34c23d027 --- /dev/null +++ b/src/Zigzag_persistence/example/zigzag_filtration_example.txt @@ -0,0 +1,49 @@ +# simple example of zigzag filtration +# i: inclusion +# r: removal +# first value: filtration value +# remaining values: if inclusion, boundary of the simplex to include in ascending order, +# a number corresponds to the id of the simplex in the boundary +# the ids start at 0 and corresponds to their "arrow number" +# If removal: id of the simplex to remove + dimension +# #: comment line + +i 0. +i 0. +i 0. + +i 1. 0 1 +i 1. 0 2 +i 1. + +i 2. 1 2 +i 2. +i 2. 5 7 + +i 3. +i 3. 3 4 6 +i 3. 7 9 +i 3. 5 9 + +i 4. 8 11 12 + +r 5. 10 2 + +r 6. 13 2 +i 6. 1 7 +i 6. 3 4 6 + +i 7. 2 7 +i 7. 8 11 12 +i 7. 0 7 +i 7. 4 18 20 +i 7. 6 16 18 +i 7. 3 16 20 + +r 8. 19 2 + +r 9. 8 1 +r 9. 12 1 +i 9. 17 21 22 23 + +r 10. 27 3 \ No newline at end of file diff --git a/src/Zigzag_persistence/include/gudhi/filtered_zigzag_persistence.h b/src/Zigzag_persistence/include/gudhi/filtered_zigzag_persistence.h new file mode 100644 index 0000000000..c026da8534 --- /dev/null +++ b/src/Zigzag_persistence/include/gudhi/filtered_zigzag_persistence.h @@ -0,0 +1,583 @@ +/* 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): Clément Maria + * + * Copyright (C) 2021 Inria + * + * Modification(s): + * - 2023/05 Hannah Schreiber: Rework of the interface, reorganization and debug + * - 2023/05 Hannah Schreiber: Addition of infinite bars + * - 2024/06 Hannah Schreiber: Separation of the zigzag algorithm from the filtration value management + * - YYYY/MM Author: Description of the modification + */ + +/** + * @file filtered_zigzag_persistence.h + * @author Clément Maria, Hannah Schreiber + * @brief Contains the implementation of the @ref Gudhi::zigzag_persistence::Default_filtered_zigzag_options structure + * and the @ref Gudhi::zigzag_persistence::Filtered_zigzag_persistence_with_storage and + * @ref Gudhi::zigzag_persistence::Filtered_zigzag_persistence classes. + */ + +#ifndef FILTERED_ZIGZAG_PERSISTENCE_H_ +#define FILTERED_ZIGZAG_PERSISTENCE_H_ + +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace Gudhi { +namespace zigzag_persistence { + +/** + * @ingroup zigzag_persistence + * + * @brief Default options for @ref Filtered_zigzag_persistence_with_storage and @ref Filtered_zigzag_persistence. + */ +struct Default_filtered_zigzag_options : Default_zigzag_options { + using Cell_key = int; /**< Cell ID used in the given boundaries. */ + using Filtration_value = double; /**< Filtration value type. */ +}; + +/** + * @ingroup zigzag_persistence + * + * @brief Class computing the zigzag persistent homology of a zigzag filtration. Algorithm based on \cite zigzag. + * Even though the insertions and removals are given in a "stream-like" way, the barcode and other values are + * stored during the whole process and not removed. It is therefore suited for smaller filtrations where the clean + * ups produce a higher overhead than the memory consumption. + * @details After construction of the class, the zigzag filtration should be given in a streaming like way, i.e., + * call @ref insert_cell, @ref remove_cell or @ref apply_identity for each step of the filtration in order of + * the filtration. To retrieve the current persistence diagram at any moment of the filtration, + * use @ref get_persistence_diagram or @ref get_index_persistence_diagram. + * + * ### Minimalistic example of usage + * + * #### Includes + * ``` + * #include + * ``` + * + * #### Useful aliases + * ``` + * using Filtered_zigzag_persistence_with_storage = Gudhi::zigzag_persistence::Filtered_zigzag_persistence_with_storage<>; + * ``` + * + * #### Construction with default values + * ``` + * Filtered_zigzag_persistence_with_storage zp; + * ``` + * + * #### Input of the zigzag sequence/filtration + * ``` + * // In all cases, it is important that the operations of insertions and removals are made **in the same order** + * // as in the zigzag filtration ones wants to compute the barcode from. + * + * // A cell can be identified in the boundaries by any given numerical label, it is just important that the given + * // filtration values are monotonous (ie., either only increasing or only decreasing). + * + * //inserts vertex 2 at filtration value 0.1 -> birth at 0.1 of 0-cycle + * zp.insert_cell(2, {}, 0, 0.1); + * //inserts vertex 4 at filtration value 0.1 -> birth at 0.1 of 0-cycle + * zp.insert_cell(4, {}, 0, 0.1); + * //inserts edge 5 = (2,4) at filtration value 0.3 -> death at 0.3 -> outputs/stores (0, 0.1, 0.3) + * zp.insert_cell(5, {2, 4}, 1, 0.3); + * //inserts vertex 3 at filtration value 0.4 -> birth at 0.4 of 0-cycle + * zp.insert_cell(3, {}, 0, 0.4); + * //inserts edge 6 = (2,3) at filtration value 0.4 -> death at 0.4 of the cycle born at 0.4 -> outputs/stores nothing + * zp.insert_cell(6, {2, 3}, 1, 0.4); + * //inserts edge 9 = (3,4) at filtration value 1.2 -> birth at 1.2 of 1-cycle + * zp.insert_cell(9, {4, 3}, 1, 1.2); + * //removes edge 6 at filtration value 1.5 -> death at 1.5 -> outputs/stores (1, 1.2, 1.5) + * zp.remove_cell(6, 1.5); + * //removes edge 5 at filtration value 2.0 -> birth at 2.0 of 0-cycle + * zp.remove_cell(5, 2.0); + * ``` + * + * #### Finalizations + * ``` + * // The bars are stored within the class and where not output at all for now. + * + * //get all bars in a vector + * auto barcode = zp.get_persistence_diagram(); + * + * //do something with the vector, e.g., stream out content: + * for (auto& bar : barcode) { + * std::cout << bar << std::endl; + * } + * ``` + * + * @tparam FilteredZigzagOptions Structure following the @ref FilteredZigzagOptions concept. + * Default value: @ref Default_filtered_zigzag_options. + */ +template +class Filtered_zigzag_persistence_with_storage +{ + public: + using Options = FilteredZigzagOptions; /**< Zigzag options. */ + using Internal_key = typename Options::Internal_key; /**< Key and index type, has to be signed. */ + using Cell_key = typename Options::Cell_key; /**< Cell ID type from external inputs. */ + using Filtration_value = typename Options::Filtration_value; /**< Type for filtration values. */ + using Dimension = typename Options::Dimension; /**< Type for dimension values. */ + + /** + * @brief Persistence index interval type. + */ + using Index_interval = Gudhi::persistence_matrix::Persistence_interval; + /** + * @brief Persistence filtration interval type. + */ + using Filtration_value_interval = Gudhi::persistence_matrix::Persistence_interval; + + /** + * @brief Constructor. + * @details After construction of the class, the zigzag filtration should be given in a streaming like way, i.e., + * call @ref insert_cell, @ref remove_cell or @ref apply_identity for each step of the filtration in order of + * the filtration. To retrieve the current persistence diagram at any moment of the filtration, + * use @ref get_persistence_diagram or @ref get_index_persistence_diagram. + * + * @param preallocationSize Reserves space for @p preallocationSize cells in the internal data structure. + * This is optional and just helps skip a few reallocations. The optimal value (no reallocation, no wasted space) is + * the number of cells in the biggest complex of the filtration. + * Default value: 0. + * @param ignoreCyclesAboveDim Ignores cycles in dimension larger or equal in the final diagram. + * If -1, no cycles are ignored. Default value: -1. + */ + Filtered_zigzag_persistence_with_storage(unsigned int preallocationSize = 0, int ignoreCyclesAboveDim = -1) + : dimMax_(ignoreCyclesAboveDim), + persistenceDiagram_(), + numArrow_(-1), + previousFiltrationValue_(std::numeric_limits::infinity()), + pers_( + [&](Dimension dim, Internal_key birth, Internal_key death) { + if (dimMax_ == -1 || (dimMax_ != -1 && dim < dimMax_)) { // don't record intervals over max dim + persistenceDiagram_.emplace_back(birth, death, dim); + } + }, + preallocationSize) {} + + /** + * @brief Updates the zigzag persistence diagram after the insertion of the given cell. + * + * @tparam BoundaryRange Range type needing size, begin and end members. + * @param cellID ID representing the inserted cell. + * @param boundary Boundary of the inserted cell. The range should be composed of the IDs of all cells contained in + * the boundary (i.e. with non-zero coefficients), using the ID specified as `cellID` when the corresponding cell + * was previously inserted (recall that the cells should be inserted in order of filtration). + * @param dimension Dimension of the inserted cell. + * @param filtrationValue Filtration value associated to the cell. + * Assumed to be always larger or equal to previously used filtration values or always smaller or equal than previous + * values, ie. the changes are monotonous. + * @return Number of the operation. + */ + template > + Internal_key insert_cell(Cell_key cellID, + const BoundaryRange& boundary, + Dimension dimension, + Filtration_value filtrationValue) + { + if (dimMax_ != -1 && dimension > dimMax_) { + return apply_identity(); + } + + ++numArrow_; + + _store_filtration_value(filtrationValue); + + [[maybe_unused]] auto res = handleToKey_.try_emplace(cellID, numArrow_); + + GUDHI_CHECK(res.second, "Zigzag_persistence::insert_cell - cell already in the complex"); + + // Compute the keys of the cells of the boundary. + std::vector translatedBoundary; + translatedBoundary.reserve(dimension * 2); // boundary does not have to have `size()` + for (auto b : boundary) { + translatedBoundary.push_back(handleToKey_.at(b)); // TODO: add possibilities of coefficients + } + std::sort(translatedBoundary.begin(), translatedBoundary.end()); + + pers_.insert_cell(translatedBoundary, dimension); + + return numArrow_; + } + + /** + * @brief Updates the zigzag persistence diagram after the removal of the given cell if the cell was contained + * in the current complex (note that it will not contain cells of dimension > ignoreCyclesAboveDim if the latter was + * non negative at construction of the class). Otherwise, just increases the operation count by one. + * + * @param cellID ID representing the cell to remove. Should be the same than the one used to insert it. + * @param filtrationValue Filtration value associated to the removal. + * Assumed to be always larger or equal to previously used filtration values or always smaller or equal than previous + * values, ie. the changes are monotonous. + * @return Number of the operation. + */ + Internal_key remove_cell(Cell_key cellID, Filtration_value filtrationValue) { + auto it = handleToKey_.find(cellID); + + if (it == handleToKey_.end()) { + return apply_identity(); + } + + ++numArrow_; + + _store_filtration_value(filtrationValue); + + pers_.remove_cell(it->second); + handleToKey_.erase(it); + + return numArrow_; + } + + /** + * @brief To use when a cell is neither inserted nor removed, but the filtration moves along the identity operator + * on homology level. Useful to keep the birth/death indices aligned when insertions/removals are purposely skipped + * to avoid useless computation. + * @return Number of the operation. + */ + Internal_key apply_identity() { + ++numArrow_; + pers_.apply_identity(); + return numArrow_; + } + + /** + * @brief Returns the "index persistence diagram" of the current filtration, that is, the pairs of atomic arrow + * numbers corresponding to a birth-death pair. Does not contain points at infinity, only the cycle classes which + * already died are represented. + * + * @return Reference to the list of intervals. + */ + const std::vector& get_index_persistence_diagram() const { return persistenceDiagram_; } + + /** + * @brief Returns the filtration value \f$f(idx)\f$ associated to the index \f$idx\f$ returned + * by @ref get_index_persistence_diagram. + * + * @param idx Birth or death index + * @return Filtration_value Filtration value associated to @p idx. + */ + Filtration_value get_filtration_value_from_index(Internal_key idx) { + // lower_bound(x) returns leftmost y s.t. x <= y + auto itBirth = + std::lower_bound(filtrationValues_.begin(), + filtrationValues_.end(), + idx, + [](std::pair p, Internal_key k) { return p.first < k; }); + if (itBirth == filtrationValues_.end() || itBirth->first > idx) { + --itBirth; + } + return itBirth->second; + }; + + /** + * @brief Returns the current persistence diagram. + * + * @param shortestInterval Threshold. Every bar shorter than the given value will not be returned. Default value: 0. + * @param includeInfiniteBars If set to true, infinite bars are included in the diagram. Default value: true. + * @return A vector of pairs of filtration values representing the persistence diagram. + */ + std::vector get_persistence_diagram(Filtration_value shortestInterval = 0., + bool includeInfiniteBars = true) { + std::vector diag = _get_persistence_diagram(shortestInterval); + + if (includeInfiniteBars) { + _retrieve_infinite_bars(diag); + } + + return diag; + } + + private: + std::unordered_map handleToKey_; /**< Map from input keys to internal keys. */ + Dimension dimMax_; /**< Maximal dimension of a bar to record. */ + std::vector persistenceDiagram_; /**< Stores current closed persistence intervals. */ + Internal_key numArrow_; /**< Current arrow number. */ + Filtration_value previousFiltrationValue_; /**< Filtration value of the previous arrow. */ + /** + * @brief filtrationValues_ stores consecutive pairs (i,f) , (j,f') with f != f', + * meaning that all inserted cells with key in [i;j-1] have filtration value f, + * i is the smallest cell index whose cell has filtration value f. + */ + std::vector > filtrationValues_; + Zigzag_persistence pers_; /**< Class computing the pairs. */ + + /** + * @brief Stores the filtration value if the value is new. Assumes that the given value is either greater (or equal) + * than previous ones or smaller (or equal) than previous ones. + * + * @param filtrationValue Filtration value to store. + */ + void _store_filtration_value(Filtration_value filtrationValue) { + if (filtrationValue != previousFiltrationValue_) // check whether the filt value has changed + { + // consecutive pairs (i,f), (j,f') mean cells of index k in [i,j-1] have filtration value f + previousFiltrationValue_ = filtrationValue; + filtrationValues_.emplace_back(numArrow_, previousFiltrationValue_); + } + } + + /** + * @brief Returns the current persistence diagram without infinite bars. + * + * @param shortestInterval Intervals shorter than the given value are ignored. + * @return Vector of intervals. + */ + std::vector _get_persistence_diagram(Filtration_value shortestInterval) { + std::vector diag; + diag.reserve(persistenceDiagram_.size()); + + // std::stable_sort(filtrationValues_.begin(), filtrationValues_.end(), + // [](std::pair p1, std::pair p2) { + // return p1.first < p2.first; + // }); + + for (auto bar : persistenceDiagram_) { + Filtration_value birth = get_filtration_value_from_index(bar.birth); + Filtration_value death = get_filtration_value_from_index(bar.death); + if (birth > death) { + std::swap(birth, death); + } + + if (death - birth > shortestInterval) { + diag.emplace_back(birth, death, bar.dim); + } + } + + return diag; + } + + /** + * @brief Computes the births of the current essential cycles. + * + * @param diag Reference to vector where to store the infinite bars. + */ + void _retrieve_infinite_bars(std::vector& diag) { + auto stream_infinite_interval = [&](Dimension dim, Internal_key birthIndex) { + if (dimMax_ == -1 || (dimMax_ != -1 && dim < dimMax_)) + diag.emplace_back(get_filtration_value_from_index(birthIndex), Filtration_value_interval::inf, dim); + }; + + pers_.get_current_infinite_intervals(stream_infinite_interval); + } +}; // end class Filtered_zigzag_persistence_with_storage + +/** + * @ingroup zigzag_persistence + * + * @brief Class computing the zigzag persistent homology of a zigzag filtration. Algorithm based on \cite zigzag. + * @details After construction of the class, the zigzag filtration should be given in a streaming like way, i.e., + * call @ref insert_cell, @ref remove_cell or @ref apply_identity for each step of the filtration in order of + * the filtration. The bars of the diagram are retrieved via the given callback method every time + * a pair with non-zero length is closed. To retrieve the open/infinite bars, use @ref get_current_infinite_intervals. + * + * ### Minimalistic example of usage + * + * #### Includes + * ``` + * #include + * ``` + * + * #### Useful aliases + * ``` + * using Filtered_zigzag_persistence = Gudhi::zigzag_persistence::Filtered_zigzag_persistence<>; + * using Dimension = Filtered_zigzag_persistence::Dimension; + * using filtration_value_type = Filtered_zigzag_persistence::Filtration_value; + * ``` + * + * #### Construction with default values + * ``` + * //Filtered_zigzag_persistence(callback) with for example callback method as a anonymous lambda + * Filtered_zigzag_persistence zp([](Dimension dim, filtration_value_type birth, filtration_value_type death) { + * std::cout << "[" << dim << "] " << birth << " - " << death << std::endl; + * }); + * ``` + * + * #### Input of the zigzag sequence/filtration + * ``` + * // In all cases, it is important that the operations of insertions and removals are made **in the same order** + * // as in the zigzag filtration ones wants to compute the barcode from. + * + * // A cell can be identified in the boundaries by any given numerical label, it is just important that the given + * // filtration values are monotonous (ie., either only increasing or only decreasing). + * + * //inserts vertex 2 at filtration value 0.1 -> birth at 0.1 of 0-cycle + * zp.insert_cell(2, {}, 0, 0.1); + * //inserts vertex 4 at filtration value 0.1 -> birth at 0.1 of 0-cycle + * zp.insert_cell(4, {}, 0, 0.1); + * //inserts edge 5 = (2,4) at filtration value 0.3 -> death at 0.3 -> outputs/stores (0, 0.1, 0.3) + * zp.insert_cell(5, {2, 4}, 1, 0.3); + * //inserts vertex 3 at filtration value 0.4 -> birth at 0.4 of 0-cycle + * zp.insert_cell(3, {}, 0, 0.4); + * //inserts edge 6 = (2,3) at filtration value 0.4 -> death at 0.4 of the cycle born at 0.4 -> outputs/stores nothing + * zp.insert_cell(6, {2, 3}, 1, 0.4); + * //inserts edge 9 = (3,4) at filtration value 1.2 -> birth at 1.2 of 1-cycle + * zp.insert_cell(9, {4, 3}, 1, 1.2); + * //removes edge 6 at filtration value 1.5 -> death at 1.5 -> outputs/stores (1, 1.2, 1.5) + * zp.remove_cell(6, 1.5); + * //removes edge 5 at filtration value 2.0 -> birth at 2.0 of 0-cycle + * zp.remove_cell(5, 2.0); + * ``` + * + * #### Finalizations + * ``` + * // Only the closed bars where output so far, so the open/infinite bars still need to be retrieved. + * + * //in this example, outputs (0, 0.1) and (0, 2.0) + * zp.get_current_infinite_intervals([](Dimension dim, filtration_value_type birth){ + * std::cout << "[" << dim << "] " << birth << " - inf" << std::endl; + * }); + * ``` + * + * @tparam FilteredZigzagOptions Structure following the @ref FilteredZigzagOptions concept. + * Default value: @ref Default_filtered_zigzag_options. + */ +template +class Filtered_zigzag_persistence { + public: + using Options = FilteredZigzagOptions; /**< Zigzag options. */ + using Internal_key = typename Options::Internal_key; /**< Key and index type, has to be signed. */ + using Cell_key = typename Options::Cell_key; /**< Cell ID type from external inputs. */ + using Filtration_value = typename Options::Filtration_value; /**< Type for filtration values. */ + using Dimension = typename Options::Dimension; /**< Type for dimension values. */ + + /** + * @brief Constructor. + * @details After construction of the class, the zigzag filtration should be given in a streaming like way, i.e., + * call @ref insert_cell, @ref remove_cell or @ref apply_identity for each step of the filtration in order of + * the filtration. The bars of the diagram are retrieved via the given callback method every time + * a pair with non-zero length is closed. To retrieve the open/infinite bars, use @ref get_current_infinite_intervals. + * + * @param stream_interval Callback method to process the birth and death values of a persistence bar. + * Has to take three arguments as input: first the dimension of the cycle, then the birth value of the cycle + * and third the death value of the cycle. The values corresponds to the filtration values which were given at + * insertions or removals. Note that bars of length 0 will not be token into account. + * @param preallocationSize Reserves space for @p preallocationSize cells in the internal data structure. + * This is optional and just helps skip a few reallocations. The optimal value (no reallocation, no wasted space) is + * the number of cells in the biggest complex of the filtration. + * Default value: 0. + * @tparam F Type of callback method. + */ + template + Filtered_zigzag_persistence(F&& stream_interval, unsigned int preallocationSize = 0) + : handleToKey_(preallocationSize), + numArrow_(-1), + keyToFiltrationValue_(preallocationSize), + pers_( + [&, stream_interval = std::forward(stream_interval)](Dimension dim, + Internal_key birth, + Internal_key death) { + auto itB = keyToFiltrationValue_.find(birth); + auto itD = keyToFiltrationValue_.find(death); + if (itB->second != itD->second) stream_interval(dim, itB->second, itD->second); + keyToFiltrationValue_.erase(itB); + keyToFiltrationValue_.erase(itD); + }, + preallocationSize) {} + + /** + * @brief Updates the zigzag persistence diagram after the insertion of the given cell. + * + * @tparam BoundaryRange Range type needing size, begin and end members. + * @param cellID ID representing the inserted cell. + * @param boundary Boundary of the inserted cell. The range should be composed of the IDs of all cells contained in + * the boundary (i.e. with non-zero coefficients), using the ID specified as `cellID` when the corresponding cell + * was previously inserted (recall that the cells should be inserted in order of filtration). + * @param dimension Dimension of the inserted cell. + * @param filtrationValue Filtration value associated to the cell. + * Assumed to be always larger or equal to previously used filtration values or always smaller or equal than previous + * values, ie. the changes are monotonous. + */ + template > + Internal_key insert_cell(Cell_key cellID, + const BoundaryRange& boundary, + Dimension dimension, + Filtration_value filtrationValue) + { + ++numArrow_; + + [[maybe_unused]] auto res = handleToKey_.try_emplace(cellID, numArrow_); + + GUDHI_CHECK(res.second, "Zigzag_persistence::insert_cell - cell already in the complex"); + + keyToFiltrationValue_.try_emplace(numArrow_, filtrationValue); + + // Compute the keys of the cells of the boundary. + std::vector translatedBoundary; + translatedBoundary.reserve(dimension * 2); // boundary does not have to have `size()` + for (auto b : boundary) { + translatedBoundary.push_back(handleToKey_.at(b)); // TODO: add possibilities of coefficients + } + std::sort(translatedBoundary.begin(), translatedBoundary.end()); + + pers_.insert_cell(translatedBoundary, dimension); + + return numArrow_; + } + + /** + * @brief Updates the zigzag persistence diagram after the removal of the given cell. + *preallocationSize + * @param cellID ID representing the cell to remove. Should be the same than the one used to insert it. + * @param filtrationValue Filtration value associated to the removal. + * Assumed to be always larger or equal to previously used filtration values or always smaller or equal than previous + * values, ie. the changes are monotonous. + */ + Internal_key remove_cell(Cell_key cellID, Filtration_value filtrationValue) { + ++numArrow_; + + auto it = handleToKey_.find(cellID); + GUDHI_CHECK(it != handleToKey_.end(), "Zigzag_persistence::remove_cell - cell not in the complex"); + + keyToFiltrationValue_.try_emplace(numArrow_, filtrationValue); + + pers_.remove_cell(it->second); + handleToKey_.erase(it); + + return numArrow_; + } + + /** + * @brief To use when a cell is neither inserted nor removed, but the filtration moves along the identity operator + * on homology level. Useful to keep the birth/death indices aligned when insertions/removals are purposely skipped + * to avoid useless computation. + */ + Internal_key apply_identity() { + ++numArrow_; + pers_.apply_identity(); + return numArrow_; + } + + /** + * @brief Outputs through the given callback method all current infinite bars. + * + * @tparam F Type of the callback method. Takes two arguments: the dimension of the cycle and the birth value + * of the cycle. + * @param stream_infinite_interval Method processing the unpaired birth values. + */ + template + void get_current_infinite_intervals(F&& stream_infinite_interval) { + pers_.get_current_infinite_intervals( + [&](Dimension dim, Internal_key birth) { stream_infinite_interval(dim, keyToFiltrationValue_.at(birth)); }); + } + + private: + template + using Dictionary = std::unordered_map; // TODO: benchmark with other map types + + Dictionary handleToKey_; /**< Map from input keys to internal keys. */ + Internal_key numArrow_; /**< Current arrow number. */ + Dictionary keyToFiltrationValue_; /**< Cell Key to filtration value map. */ + Zigzag_persistence pers_; /**< Class computing the pairs. */ +}; // end class Filtered_zigzag_persistence + +} // namespace zigzag_persistence +} // namespace Gudhi + +#endif // FILTERED_ZIGZAG_PERSISTENCE_H_ diff --git a/src/Zigzag_persistence/include/gudhi/zigzag_persistence.h b/src/Zigzag_persistence/include/gudhi/zigzag_persistence.h new file mode 100644 index 0000000000..e704587dee --- /dev/null +++ b/src/Zigzag_persistence/include/gudhi/zigzag_persistence.h @@ -0,0 +1,488 @@ +/* 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): Clément Maria + * + * Copyright (C) 2021 Inria + * + * Modification(s): + * - 2023/05 Hannah Schreiber: Rework of the interface, reorganization and debug + * - 2023/05 Hannah Schreiber: Addition of infinite bars + * - 2024/06 Hannah Schreiber: Separation of the zigzag algorithm from the filtration value management + * - YYYY/MM Author: Description of the modification + */ + +/** + * @file zigzag_persistence.h + * @author Clément Maria, Hannah Schreiber + * @brief Contains the implementation of the @ref Gudhi::zigzag_persistence::Zigzag_persistence class. + */ + +#ifndef ZIGZAG_PERSISTENCE_H_ +#define ZIGZAG_PERSISTENCE_H_ + +#include +#include +#include +#include + +#include +#if BOOST_VERSION >= 108100 +#include //don't exist for lower versions of boost +// #include +#else +#include +#endif + +#include + +namespace Gudhi { +namespace zigzag_persistence { + +/** + * @ingroup zigzag_persistence + * + * @brief Options for the internal matrix of @ref Zigzag_persistence. + * + * @tparam column_type Column type of the matrix. + */ +template +struct Zigzag_matrix_options : Gudhi::persistence_matrix::Default_options { + static const bool has_row_access = true; + static const bool has_column_pairings = false; + static const bool has_vine_update = true; + static const bool is_of_boundary_type = false; + static const bool has_map_column_container = true; + static const bool has_removable_columns = true; + static const bool has_removable_rows = true; +}; + +/** + * @ingroup zigzag_persistence + * + * @brief Default options for @ref Zigzag_persistence. + */ +struct Default_zigzag_options { + using Internal_key = int; /**< Cell ID used internally, must be signed. */ + using Dimension = int; /**< Dimension value type. */ + /** + * @brief Column type use by the internal matrix. + */ + static const Gudhi::persistence_matrix::Column_types column_type = + Gudhi::persistence_matrix::Column_types::NAIVE_VECTOR; +}; + +// TODO: add the possibility of something else than Z2. Which means that the possibility of vineyards without Z2 +// also needs to be implemented. The theory needs to be done first. +/** + * @class Zigzag_persistence zigzag_persistence.h gudhi/zigzag_persistence.h + * @brief Class computing the zigzag persistent homology of a zigzag sequence. Algorithm based on \cite zigzag. + * @details After construction of the class, the zigzag filtration should be given in a streaming like way, i.e., + * call @ref insert_cell, @ref remove_cell or @ref apply_identity for each step of the filtration in order of + * the filtration. The pairs of birth and death indices are retrieved via the given callback method every time + * a pair is closed. To retrieve the open pairs (corresponding to infinite bars), + * use @ref get_current_infinite_intervals. + * + * ### Minimalistic example of usage + * + * #### Includes + * ``` + * #include + * ``` + * + * #### Useful aliases + * ``` + * using Zigzag_persistence = Gudhi::zigzag_persistence::Zigzag_persistence<>; + * + * using Dimension = Zigzag_persistence::Dimension; + * using Index = Zigzag_persistence::Index; + * ``` + * + * #### Construction with default values + * ``` + * //Zigzag_persistence(callback) with for example callback method as a anonymous lambda + * Zigzag_persistence zp([](Dimension dim, Index birth, Index death) { + * std::cout << "[" << dim << "] " << birth << " - " << death << std::endl; + * }); + * ``` + * + * #### Input of the zigzag sequence/filtration + * ``` + * // In all cases, it is important that the operations of insertions and removals are made **in the same order** + * // as in the zigzag filtration ones wants to compute the barcode from. + * + * // A cell has to be identified in the boundaries by the operation number the cell was inserted with in the sequence. + * + * //inserts vertex 0 -> birth at 0 of 0-cycle + * zp.insert_cell({}, 0); + * //inserts vertex 1 -> birth at 1 of 0-cycle + * zp.insert_cell({}, 0); + * //inserts edge 2 = (0,1) -> death at 2 -> outputs (0, 1, 2) + * zp.insert_cell({0, 1}, 1); + * //inserts vertex 3 -> birth at 3 of 0-cycle + * zp.insert_cell({}, 0); + * //inserts edge 4 = (0,3) -> death at 4 -> outputs (0, 3, 4) + * zp.insert_cell({0, 3}, 1); + * //inserts edge 5 = (1,3) -> birth at 5 of 1-cycle + * zp.insert_cell({1, 3}, 1); + * //removes edge 4 -> death at 6 -> outputs (1, 5, 6) + * zp.remove_cell(4); + * //removes edge 2 -> birth at 7 of 0-cycle + * zp.remove_cell(2); + * ``` + * + * #### Finalizations + * ``` + * // Only the closed bars where output so far, so the open/infinite bars still need to be retrieved. + * + * //in this example, outputs (0, 0) and (0, 7) + * zp.get_current_infinite_intervals([](Dimension dim, Index birth){ + * std::cout << "[" << dim << "] " << birth << " - inf" << std::endl; + * }); + * ``` + * + * @ingroup zigzag_persistence + * + * @tparam ZigzagOptions Structure following the @ref ZigzagOptions concept. Default value: @ref Default_zigzag_options. + */ +template +class Zigzag_persistence +{ + public: + using Options = ZigzagOptions; /**< Zigzag options. */ + using Index = typename Options::Internal_key; /**< Key and index type, has to be signed. */ + using Dimension = typename Options::Dimension; /**< Type for dimension values. */ + + private: +#if BOOST_VERSION >= 108100 + using Birth_dictionary = boost::unordered_flat_map; /**< Dictionary type. */ + // using Birth_dictionary = boost::unordered_map; /**< Dictionary type. */ +#else + using Birth_dictionary = std::unordered_map; /**< Dictionary type. */ +#endif + using Matrix_options = Zigzag_matrix_options; /**< Matrix options. */ + using Matrix = Gudhi::persistence_matrix::Matrix; /**< Matrix. */ + using Matrix_index = typename Matrix::Index; /**< Matrix indexation type. */ + + /** \brief Maintains the birth ordering \f$\leq_b\f$. + * + * \details Contains a map of size the number of + * non-zero rows of the homology matrix, at any time during the computation of + * zigzag persistence. + * + * By construction, we maintain the map satisfying + * 'birthToPos_[i] < birthToPos_[j]', with \f$0 <= i,j <= k\f$ indices in the quiver + * '\f$0 \leftrightarrow ... \leftrightarrow i \leftrightarrow ... \leftrightarrow k\f$' visited at time + * \f$k\f$ of the algorithm (prefix of length \f$k\f$ of the full zigzag filtration + * '\f$0 \leftrightarrow ... \leftrightarrow i \leftrightarrow ... + * \leftrightarrow k \leftrightarrow ... \leftrightarrow n\f$' + * that is studied), iff \f$i <_b j\f$ for the birth ordering. + * + * By construction, when adding index \f$k+1\f$ to '\f$0 \leftrightarrow ... \leftrightarrow i \leftrightarrow .. + * \leftrightarrow k \leftrightarrow k+1'\f$, we have: + * - if \f$k \rightarrow k+1\f$ forward, then \f$j <_b k+1\f$ for all indices \f$j < k+1\f$, otherwise + * - if \f$k \leftarrow k+1\f$ backward, then \f$k+1 <_b j\f$ for all indices \f$j < k+1\f$. + */ + struct Birth_ordering { + /** + * @brief Default constructor + */ + Birth_ordering() : birthToPos_(), maxBirthPos_(0), minBirthPos_(-1) {} + + /** + * @brief Inserts arrow number in the ordering after an insertion. + * When the arrow key-1 -> key is forward, key is larger than any other index + * i < key in the birth ordering b k2. + * + * @param k1 + * @param k2 + * @return true if k1 >b k2, false otherwise. + */ + bool reverse_birth_order(Index k1, Index k2) const { return birthToPos_.at(k1) > birthToPos_.at(k2); } + + private: + Birth_dictionary birthToPos_; /**< birth_to_pos_[i] < birth_to_pos_[j] iff i stream_interval, + unsigned int preallocationSize = 0) + : matrix_( + preallocationSize, + [this](Matrix_index columnIndex1, Matrix_index columnIndex2) -> bool { + if (matrix_.get_column(columnIndex1).is_paired()) { + return matrix_.get_pivot(columnIndex1) < matrix_.get_pivot(columnIndex2); + } + return birthOrdering_.birth_order(births_.at(columnIndex1), births_.at(columnIndex2)); + }, + [this](Matrix_index columnIndex1, Matrix_index columnIndex2) -> bool { return false; }), + numArrow_(-1), + stream_interval_(std::move(stream_interval)) {} + /** + * @brief Updates the zigzag persistence diagram after the insertion of the given cell. + * + * @tparam BoundaryRange Range type needing size, begin and end members. + * @param boundary Boundary of the inserted cell. The boundary should be represented by all the cells with + * non-zero coefficients generating it. A cell should be represented by the arrow number when the cell appeared for + * the first time in the filtration (if a cell was inserted and then removed and reinserted etc., only the last + * insertion counts). The cell range should be ordered by increasing arrow numbers. + * @param dimension Dimension of the inserted cell. + * @return Number of the operation. + */ + template > + Index insert_cell(const BoundaryRange& boundary, Dimension dimension) { + ++numArrow_; + _process_forward_arrow(boundary, dimension); + return numArrow_; + } + + /** + * @brief Updates the zigzag persistence diagram after the removal of the given cell. + * + * @param arrowNumber Arrow number of when the cell to remove was inserted for the last time. + * @return Number of the operation. + */ + Index remove_cell(Index arrowNumber) { + ++numArrow_; + _process_backward_arrow(arrowNumber); + return numArrow_; + } + + /** + * @brief To use when a cell is neither inserted nor removed, but the filtration moves along the identity operator + * on homology level. Useful to keep the birth/death indices aligned when insertions/removals are purposely skipped + * to avoid useless computation. Increases the arrow number by one. + * @return Number of the operation. + */ + Index apply_identity() { return ++numArrow_; } + + /** + * @brief Outputs through the given callback method all birth indices which are currently not paired with + * a death index. + * + * @tparam F Type of the callback method. Takes two arguments: the dimension of the cycle and the birth index + * of the cycle. + * @param stream_infinite_interval Method processing the unpaired birth indices. + */ + template + void get_current_infinite_intervals(F&& stream_infinite_interval) { + for (auto& p : births_) { + auto& col = matrix_.get_column(p.first); + stream_infinite_interval(col.get_dimension(), p.second); + } + } + + private: + /** + * @brief Express the boundary cycle of the new cell as a sum of cycles in a matrix. + * If some cycles are not boundary cycles, i.e., columns with F-index + * in the matrix, it applies a surjective diamond to the zigzag module. + * + * @param boundary Boundary of the inserted cell. + * @param dim Dimension of the inserted cell. + */ + template + void _process_forward_arrow(const BoundaryRange& boundary, Dimension dim) { + std::vector chainsInF = matrix_.insert_boundary(numArrow_, boundary, dim); + + if (!chainsInF.empty()) { + _apply_surjective_reflection_diamond(dim, chainsInF); + } else { + birthOrdering_.add_birth_forward(numArrow_); + births_[matrix_.get_column_with_pivot(numArrow_)] = numArrow_; + } + } + + /** + * @brief Applies the surjective reflection diamond principle to the current filtration. + * + * @details The vector chainsInF is sorted by decreasing lowest index values in the + * columns corresponding to the chains, due to its computation in the reduction of + * the boundary in _process_forward_arrow(...). It is equivalent to decreasing death index + * order w.r.t. the & chainsInF) { + // fp is the largest death index for <=d + // Set col_fp: col_fp <- col_f1+...+col_fp (now in G); preserves lowest idx + auto chainFp = chainsInF[0]; // col_fp, with largest death bool { + return birthOrdering_.reverse_birth_order(k1, k2); + }; // true iff b(k1) >b b(k2) + + // availableBirth: for all i by >d value of the d_i, + // contains at step i all b_j, j > i, and maybe b_i if not stolen + std::set availableBirth(cmp_birth); + // for f1 to f_{p} (i by <=d), insertion in availableBirth sorts by >=b + for (auto& chainF : chainsInF) { + availableBirth.insert(births_.at(chainF)); + } + + auto maxbIt = availableBirth.begin(); // max birth cycle + auto maxb = *maxbIt; // max birth value, for persistence diagram + availableBirth.erase(maxbIt); // remove max birth cycle (stolen) + + auto lastModifiedChainIt = chainsInF.rbegin(); + + // consider all death indices by increasing the maximal available death. + // Let c_1 ... c_f be the chains s.t. <[c_1+...+c_f]> is the kernel and + // death(c_i) >d death(c_i-1). If the birth of c_i is not available, we set + // c_i <- c_i + c_i-1 + ... + c_1, which is [c_i + c_i-1 + ... + c_1] on + // the right (of death the maximali <=> the maxj>k, are indeed c_j + // set c_i <- c_i + (c_i-1) + ... + (c_k+1) + (c_k + ... + c_1) + for (auto chainPassedIt = lastModifiedChainIt; chainPassedIt != chainFIt; ++chainPassedIt) { + // all with smaller modifiedColumns; + const auto& row = matrix_.get_row(cellID); + modifiedColumns.reserve(row.size()); + std::transform(row.begin(), row.end(), std::back_inserter(modifiedColumns), + [](const auto& cell) { return cell.get_column_index(); }); + // Sort by left-to-right order in the matrix_ (no order maintained in rows) + std::stable_sort(modifiedColumns.begin(), modifiedColumns.end(), [this](Matrix_index i1, Matrix_index i2) { + return matrix_.get_pivot(i1) < matrix_.get_pivot(i2); + }); + + // Modifies curr_col, not the other one. + for (auto otherColIt = std::next(modifiedColumns.begin()); otherColIt != modifiedColumns.end(); ++otherColIt) { + currCol = matrix_.vine_swap_with_z_eq_1_case(currCol, *otherColIt); + } + + // curr_col points to the column to remove by restriction of K to K-{\sigma} + auto& col = matrix_.get_column(currCol); + if (!col.is_paired()) { // in F + auto it = births_.find(currCol); + stream_interval_(col.get_dimension(), it->second, numArrow_); + birthOrdering_.remove_birth(it->second); + births_.erase(it); + } else { // in H -> paired with c_g, that now belongs to F now + // maintain the <=b order + birthOrdering_.add_birth_backward(numArrow_); + births_[col.get_paired_chain_index()] = numArrow_; + } + + // cannot be in G as the removed cell is maximal + matrix_.remove_maximal_cell(cellID, {}); // also un-pairs c_g if in H + } + + private: + Matrix matrix_; /**< Matrix storing a base of the current chain complex. */ + Birth_dictionary births_; /**< Map cell index in F to corresponding birth. */ + Birth_ordering birthOrdering_; /**< Maintains stream_interval_; /**< Callback method for closed pairs. */ +}; // end class Zigzag_persistence + +} // namespace zigzag_persistence +} // namespace Gudhi + +#endif // ZIGZAG_PERSISTENCE_H_ diff --git a/src/Zigzag_persistence/test/CMakeLists.txt b/src/Zigzag_persistence/test/CMakeLists.txt new file mode 100644 index 0000000000..ee0fe2c4e9 --- /dev/null +++ b/src/Zigzag_persistence/test/CMakeLists.txt @@ -0,0 +1,12 @@ +project(Zigzag_persistence_tests) + +include(GUDHI_boost_test) + +add_executable_with_targets(Zigzag_persistence_test_zigzag_persistence test_zigzag_persistence.cpp TBB::tbb) +gudhi_add_boost_test(Zigzag_persistence_test_zigzag_persistence) + +add_executable_with_targets(Zigzag_persistence_test_filtered_zigzag_persistence test_filtered_zigzag_persistence.cpp TBB::tbb) +gudhi_add_boost_test(Zigzag_persistence_test_filtered_zigzag_persistence) + + + diff --git a/src/Zigzag_persistence/test/test_filtered_zigzag_persistence.cpp b/src/Zigzag_persistence/test/test_filtered_zigzag_persistence.cpp new file mode 100644 index 0000000000..6d03c0cc9e --- /dev/null +++ b/src/Zigzag_persistence/test/test_filtered_zigzag_persistence.cpp @@ -0,0 +1,473 @@ +/* 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): Hannah Schreiber + * + * Copyright (C) 2023 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include +#include +#include + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "zigzag_persistence" +#include + +#include + +struct Interval_comparator { + Interval_comparator() {} + template + bool operator()(const Interval_filtration& p, const Interval_filtration& q) { + if (p.dim != q.dim) { + return p.dim < q.dim; + } + if (p.birth != q.birth) { + return p.birth < q.birth; + } + return p.death < q.death; + } +}; + +BOOST_AUTO_TEST_CASE(constructor) { + using ZP1 = Gudhi::zigzag_persistence::Filtered_zigzag_persistence_with_storage<>; + BOOST_CHECK_NO_THROW(ZP1 zp); + BOOST_CHECK_NO_THROW(ZP1 zp(28)); + BOOST_CHECK_NO_THROW(ZP1 zp(28, 2)); + ZP1 zp1; + BOOST_CHECK(zp1.get_persistence_diagram(0).empty()); + + using ZP2 = Gudhi::zigzag_persistence::Filtered_zigzag_persistence<>; + std::vector > pairs; + auto stream = [&](int dim, double birth, double death){ pairs.emplace_back(dim, birth, death); }; + BOOST_CHECK_NO_THROW(ZP2 zp(stream)); + BOOST_CHECK_NO_THROW(ZP2 zp(stream, 28)); + + ZP2 zp2(stream); + BOOST_CHECK(pairs.empty()); +} + +template +void test_barcode(ZP& zp, std::vector& barcode) { + auto bars = zp.get_persistence_diagram(0, true); + std::stable_sort(bars.begin(), bars.end(), Interval_comparator()); + std::stable_sort(barcode.begin(), barcode.end(), Interval_comparator()); + auto it = barcode.begin(); + for (const auto& interval : bars) { + BOOST_CHECK_EQUAL(interval.dim, it->dim); + BOOST_CHECK_EQUAL(interval.birth, it->birth); + BOOST_CHECK_EQUAL(interval.death, it->death); + ++it; + } + BOOST_CHECK(it == barcode.end()); +} + +template +void test_indices(ZP& zp, std::vector& indices, + std::vector& indexToFil) { + auto it = indices.begin(); + for (const auto& interval : zp.get_index_persistence_diagram()) { + BOOST_CHECK_EQUAL(interval.dim, it->dim); + BOOST_CHECK_EQUAL(interval.birth, it->birth); + BOOST_CHECK_EQUAL(interval.death, it->death); + BOOST_CHECK_EQUAL(zp.get_filtration_value_from_index(interval.birth), indexToFil[interval.birth]); + BOOST_CHECK_EQUAL(zp.get_filtration_value_from_index(interval.death), indexToFil[interval.death]); + ++it; + } + BOOST_CHECK(it == indices.end()); +} + +std::vector > get_boundaries() { + return {{}, + {}, + {}, + {0, 1}, + {0, 2}, + {}, + {1, 2}, + {}, + {5, 7}, + {}, + {3, 4, 6}, + {7, 9}, + {5, 9}, + {8, 11, 12}, + {10}, // remove + {13}, // remove + {1, 7}, + {3, 4, 6}, + {2, 7}, + {8, 11, 12}, + {0, 7}, + {4, 18, 20}, + {6, 16, 18}, + {3, 16, 20}, + {19}, // remove + {8}, // remove + {12}, // remove + {17, 21, 22, 23}, + {27}}; // remove +} + +std::vector get_filtration_values() { + return {0, 0, 0, + 1, 1, 1, + 2, 2, 2, + 3, 3, 3, 3, + 4, + 5, + 6, 6, 6, + 7, 7, 7, 7, 7, 7, + 8, + 9, 9, 9, + 10}; +} + +template +void test_filtered_zigzag_with_storage() { + using cell_handle = typename ZP::Cell_key; + using Filtration_value = typename ZP::Filtration_value; + using Interval_index = typename ZP::Index_interval; + using Interval_filtration = typename ZP::Filtration_value_interval; + + ZP zp(28); + std::vector realIndices; + std::vector realBarcode; + realIndices.reserve(13); + realBarcode.reserve(9); + + std::vector > simplices = get_boundaries(); + std::vector filValues = get_filtration_values(); + + for (unsigned int i = 0; i < 14; ++i) { + zp.insert_cell(i, simplices[i], simplices[i].size() == 0 ? 0 : simplices[i].size() - 1, filValues[i]); + } + + realIndices.emplace_back(1, 3, 0); + realIndices.emplace_back(2, 4, 0); + realIndices.emplace_back(7, 8, 0); + realIndices.emplace_back(6, 10, 1); + realIndices.emplace_back(9, 11, 0); + realIndices.emplace_back(12, 13, 1); + + realBarcode.emplace_back(0, 1, 0); + realBarcode.emplace_back(0, 1, 0); + realBarcode.emplace_back(2, 3, 1); + realBarcode.emplace_back(3, 4, 1); + + for (unsigned int i = 14; i < 16; ++i) { + auto id = simplices[i][0]; + zp.remove_cell(id, filValues[i]); + } + + for (unsigned int i = 16; i < 24; ++i) { + zp.insert_cell(i, simplices[i], simplices[i].size() == 0 ? 0 : simplices[i].size() - 1, filValues[i]); + } + + realIndices.emplace_back(5, 16, 0); + realIndices.emplace_back(14, 17, 1); + realIndices.emplace_back(15, 19, 1); + realIndices.emplace_back(20, 21, 1); + realIndices.emplace_back(18, 22, 1); + + realBarcode.emplace_back(1, 6, 0); + realBarcode.emplace_back(5, 6, 1); + realBarcode.emplace_back(6, 7, 1); + + for (unsigned int i = 24; i < 27; ++i) { + auto id = simplices[i][0]; + zp.remove_cell(id, filValues[i]); + } + + realIndices.emplace_back(24, 25, 1); + realBarcode.emplace_back(8, 9, 1); + + zp.insert_cell(27, simplices[27], simplices[27].size() == 0 ? 0 : simplices[27].size() - 1, filValues[27]); + + realIndices.emplace_back(23, 27, 2); + realBarcode.emplace_back(7, 9, 2); + + auto id = simplices[28][0]; + zp.remove_cell(id, filValues[28]); + + realBarcode.emplace_back(0, Interval_filtration::inf, 0); + realBarcode.emplace_back(9, Interval_filtration::inf, 0); + realBarcode.emplace_back(10, Interval_filtration::inf, 2); + + test_indices(zp, realIndices, filValues); + test_barcode(zp, realBarcode); +} + +template +void test_filtered_zigzag_with_storage_max1() { + using cell_handle = typename ZP::Cell_key; + using Filtration_value = typename ZP::Filtration_value; + using Interval_index = typename ZP::Index_interval; + using Interval_filtration = typename ZP::Filtration_value_interval; + + ZP zp(28, 1); + std::vector realIndices; + std::vector realBarcode; + realIndices.reserve(5); + realBarcode.reserve(3); + + std::vector > simplices = get_boundaries(); + std::vector filValues = get_filtration_values(); + + for (unsigned int i = 0; i < 14; ++i) { + zp.insert_cell(i, simplices[i], simplices[i].size() == 0 ? 0 : simplices[i].size() - 1, filValues[i]); + } + + realIndices.emplace_back(1, 3, 0); + realIndices.emplace_back(2, 4, 0); + realIndices.emplace_back(7, 8, 0); + realIndices.emplace_back(9, 11, 0); + + realBarcode.emplace_back( 0, 1, 0); + realBarcode.emplace_back(0, 1, 0); + + for (unsigned int i = 14; i < 16; ++i) { + auto id = simplices[i][0]; + zp.remove_cell(id, filValues[i]); + } + + for (unsigned int i = 16; i < 24; ++i) { + zp.insert_cell(i, simplices[i], simplices[i].size() == 0 ? 0 : simplices[i].size() - 1, filValues[i]); + } + + realIndices.emplace_back(5, 16, 0); + realBarcode.emplace_back(1, 6, 0); + + for (unsigned int i = 24; i < 27; ++i) { + auto id = simplices[i][0]; + zp.remove_cell(id, filValues[i]); + } + + zp.insert_cell(27, simplices[27], simplices[27].size() == 0 ? 0 : simplices[27].size() - 1, filValues[27]); + auto id = simplices[28][0]; + zp.remove_cell(id, filValues[28]); + + realBarcode.emplace_back(0, Interval_filtration::inf, 0); + realBarcode.emplace_back(9, Interval_filtration::inf, 0); + + test_indices(zp, realIndices, filValues); + test_barcode(zp, realBarcode); +} + +BOOST_AUTO_TEST_CASE(filtered_zigzag_persistence_with_storage) { + test_filtered_zigzag_with_storage >(); + test_filtered_zigzag_with_storage_max1 >(); +} + +template +void test_filtered_zigzag() { + using cell_handle = typename ZP::Cell_key; + using Filtration_value = typename ZP::Filtration_value; + using Dimension = typename ZP::Dimension; + using Interval = std::tuple; + + Interval interval; + ZP zp([&](Dimension dim, Filtration_value birth, Filtration_value death){ + BOOST_CHECK_EQUAL(std::get<0>(interval), dim); + BOOST_CHECK_EQUAL(std::get<1>(interval), birth); + BOOST_CHECK_EQUAL(std::get<2>(interval), death); + },28); + + std::vector realBarcode; + realBarcode.reserve(28); + realBarcode.emplace_back(3, 0, 0); //dummy + realBarcode.emplace_back(3, 0, 1); //dummy + realBarcode.emplace_back(3, 0, 2); //dummy + realBarcode.emplace_back(0, 0, 1); //1-3 + realBarcode.emplace_back(0, 0, 1); //2-4 + realBarcode.emplace_back(3, 0, 5); //dummy + realBarcode.emplace_back(3, 0, 6); //dummy + realBarcode.emplace_back(3, 0, 7); //dummy + realBarcode.emplace_back(0, 7, 8); //dummy + realBarcode.emplace_back(3, 0, 9); //dummy + realBarcode.emplace_back(1, 2, 3); //6-10 + realBarcode.emplace_back(0, 9, 11); //dummy + realBarcode.emplace_back(3, 0, 12); //dummy + realBarcode.emplace_back(1, 3, 4); //12-13 + realBarcode.emplace_back(3, 0, 14); //dummy + realBarcode.emplace_back(3, 0, 15); //dummy + realBarcode.emplace_back(0, 1, 6); //5-16 + realBarcode.emplace_back(1, 5, 6); //14-17 + realBarcode.emplace_back(3, 0, 18); //dummy + realBarcode.emplace_back(1, 6, 7); //15-19 + realBarcode.emplace_back(3, 0, 20); //dummy + realBarcode.emplace_back(1, 20, 21); //dummy + realBarcode.emplace_back(1, 18, 22); //dummy + realBarcode.emplace_back(3, 0, 23); //dummy + realBarcode.emplace_back(3, 0, 24); //dummy + realBarcode.emplace_back(1, 8, 9); //24-25 + realBarcode.emplace_back(3, 0, 26); //dummy + realBarcode.emplace_back(2, 7, 9); //23-27 + realBarcode.emplace_back(3, 0, 28); //dummy + + std::vector > simplices = get_boundaries(); + std::vector filValues = get_filtration_values(); + + for (unsigned int i = 0; i < 14; ++i) { + interval = realBarcode[i]; + zp.insert_cell(i, simplices[i], simplices[i].size() == 0 ? 0 : simplices[i].size() - 1, filValues[i]); + } + + for (unsigned int i = 14; i < 16; ++i) { + interval = realBarcode[i]; + auto id = simplices[i][0]; + zp.remove_cell(id, filValues[i]); + } + + for (unsigned int i = 16; i < 24; ++i) { + interval = realBarcode[i]; + zp.insert_cell(i, simplices[i], simplices[i].size() == 0 ? 0 : simplices[i].size() - 1, filValues[i]); + } + + for (unsigned int i = 24; i < 27; ++i) { + interval = realBarcode[i]; + auto id = simplices[i][0]; + zp.remove_cell(id, filValues[i]); + } + + interval = realBarcode[27]; + zp.insert_cell(27, simplices[27], simplices[27].size() == 0 ? 0 : simplices[27].size() - 1, filValues[27]); + + interval = realBarcode[28]; + auto id = simplices[28][0]; + zp.remove_cell(id, filValues[28]); + + //there is no real guarantee on the order of the infinite bars + std::vector infiniteBars; + zp.get_current_infinite_intervals([&](Dimension dim, Filtration_value birth) { + infiniteBars.emplace_back(dim, birth, std::numeric_limits::infinity()); + }); + + realBarcode.clear(); + realBarcode.emplace_back(0, 0, std::numeric_limits::infinity()); + realBarcode.emplace_back(0, 9, std::numeric_limits::infinity()); + realBarcode.emplace_back(2, 10, std::numeric_limits::infinity()); + + std::sort(infiniteBars.begin(), infiniteBars.end()); + std::sort(realBarcode.begin(), realBarcode.end()); + auto it = realBarcode.begin(); + for (const auto& interval : infiniteBars) { + BOOST_CHECK_EQUAL(std::get<0>(interval), std::get<0>(*it)); + BOOST_CHECK_EQUAL(std::get<1>(interval), std::get<1>(*it)); + ++it; + } + BOOST_CHECK(it == realBarcode.end()); +} + +template +void test_filtered_zigzag_max1() { + using cell_handle = typename ZP::Cell_key; + using Filtration_value = typename ZP::Filtration_value; + using Dimension = typename ZP::Dimension; + using Interval = std::tuple; + + Interval interval; + ZP zp([&](Dimension dim, Filtration_value birth, Filtration_value death){ + if (dim < 1){ + BOOST_CHECK_EQUAL(std::get<0>(interval), dim); + BOOST_CHECK_EQUAL(std::get<1>(interval), birth); + BOOST_CHECK_EQUAL(std::get<2>(interval), death); + } else { + BOOST_CHECK_NE(std::get<0>(interval), 0); + } + },28); + + std::vector realBarcode; + realBarcode.reserve(28); + realBarcode.emplace_back(1, 0, 0); //dummy + realBarcode.emplace_back(1, 0, 1); //dummy + realBarcode.emplace_back(1, 0, 2); //dummy + realBarcode.emplace_back(0, 0, 1); //1-3 + realBarcode.emplace_back(0, 0, 1); //2-4 + realBarcode.emplace_back(1, 0, 5); //dummy + realBarcode.emplace_back(1, 0, 6); //dummy + realBarcode.emplace_back(1, 0, 7); //dummy + realBarcode.emplace_back(1, 7, 8); //dummy + realBarcode.emplace_back(1, 0, 9); //dummy + realBarcode.emplace_back(1, 2, 3); //6-10 + realBarcode.emplace_back(1, 9, 11); //dummy + realBarcode.emplace_back(1, 0, 12); //dummy + realBarcode.emplace_back(1, 3, 4); //12-13 + realBarcode.emplace_back(1, 0, 14); //dummy + realBarcode.emplace_back(1, 0, 15); //dummy + realBarcode.emplace_back(0, 1, 6); //5-16 + realBarcode.emplace_back(1, 5, 6); //14-17 + realBarcode.emplace_back(1, 0, 18); //dummy + realBarcode.emplace_back(1, 6, 7); //15-19 + realBarcode.emplace_back(1, 0, 20); //dummy + realBarcode.emplace_back(1, 20, 21); //dummy + realBarcode.emplace_back(1, 18, 22); //dummy + realBarcode.emplace_back(1, 0, 23); //dummy + realBarcode.emplace_back(1, 0, 24); //dummy + realBarcode.emplace_back(1, 8, 9); //24-25 + realBarcode.emplace_back(1, 0, 26); //dummy + realBarcode.emplace_back(2, 7, 9); //23-27 + realBarcode.emplace_back(1, 0, 28); //dummy + + std::vector > simplices = get_boundaries(); + std::vector filValues = get_filtration_values(); + + for (unsigned int i = 0; i < 14; ++i) { + interval = realBarcode[i]; + zp.insert_cell(i, simplices[i], simplices[i].size() == 0 ? 0 : simplices[i].size() - 1, filValues[i]); + } + + for (unsigned int i = 14; i < 16; ++i) { + interval = realBarcode[i]; + auto id = simplices[i][0]; + zp.remove_cell(id, filValues[i]); + } + + for (unsigned int i = 16; i < 24; ++i) { + interval = realBarcode[i]; + zp.insert_cell(i, simplices[i], simplices[i].size() == 0 ? 0 : simplices[i].size() - 1, filValues[i]); + } + + for (unsigned int i = 24; i < 27; ++i) { + interval = realBarcode[i]; + auto id = simplices[i][0]; + zp.remove_cell(id, filValues[i]); + } + + interval = realBarcode[27]; + zp.insert_cell(27, simplices[27], simplices[27].size() == 0 ? 0 : simplices[27].size() - 1, filValues[27]); + + interval = realBarcode[28]; + auto id = simplices[28][0]; + zp.remove_cell(id, filValues[28]); + + //there is no real guarantee on the order of the infinite bars + std::vector infiniteBars; + zp.get_current_infinite_intervals([&](Dimension dim, Filtration_value birth) { + if (dim < 1){ + infiniteBars.emplace_back(dim, birth, std::numeric_limits::infinity()); + } + }); + + realBarcode.clear(); + realBarcode.emplace_back(0, 0, std::numeric_limits::infinity()); + realBarcode.emplace_back(0, 9, std::numeric_limits::infinity()); + + std::sort(infiniteBars.begin(), infiniteBars.end()); + std::sort(realBarcode.begin(), realBarcode.end()); + auto it = realBarcode.begin(); + for (const auto& interval : infiniteBars) { + BOOST_CHECK_EQUAL(std::get<0>(interval), std::get<0>(*it)); + BOOST_CHECK_EQUAL(std::get<1>(interval), std::get<1>(*it)); + ++it; + } + BOOST_CHECK(it == realBarcode.end()); +} + +BOOST_AUTO_TEST_CASE(filtered_zigzag_persistence) { + test_filtered_zigzag >(); + test_filtered_zigzag_max1 >(); +} diff --git a/src/Zigzag_persistence/test/test_zigzag_persistence.cpp b/src/Zigzag_persistence/test/test_zigzag_persistence.cpp new file mode 100644 index 0000000000..0f1edb9a2b --- /dev/null +++ b/src/Zigzag_persistence/test/test_zigzag_persistence.cpp @@ -0,0 +1,206 @@ +/* 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): Hannah Schreiber + * + * Copyright (C) 2023 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "zigzag_persistence" +#include + +#include + +using ZP = Gudhi::zigzag_persistence::Zigzag_persistence<>; + +struct Interval { + Interval() {} + Interval(int dim, ZP::Index b, ZP::Index d) : dim_(dim), b_(b), d_(d) {} + + int dim() const { return dim_; } + int birth() const { return b_; } + int death() const { return d_; } + +private: + int dim_; + ZP::Index b_; + ZP::Index d_; +}; + +BOOST_AUTO_TEST_CASE(constructor) { + std::vector pairs; + auto stream = [&](int dim, ZP::Index birth, ZP::Index death){ pairs.emplace_back(dim, birth, death); }; + BOOST_CHECK_NO_THROW(ZP zp(stream)); + BOOST_CHECK_NO_THROW(ZP zp(stream, 28)); + + ZP zp(stream); + BOOST_CHECK(pairs.empty()); +} + +void test_indices(std::vector& zp_indices, std::vector& witness_indices) { + auto it = witness_indices.begin(); + for (const Interval& interval : zp_indices) { + BOOST_CHECK_EQUAL(interval.dim(), it->dim()); + BOOST_CHECK_EQUAL(interval.birth(), it->birth()); + BOOST_CHECK_EQUAL(interval.death(), it->death()); + ++it; + } + BOOST_CHECK(it == witness_indices.end()); +} + +std::vector > get_boundaries() { + return {{}, + {}, + {}, + {0, 1}, + {0, 2}, + {}, + {1, 2}, + {}, + {5, 7}, + {}, + {3, 4, 6}, + {7, 9}, + {5, 9}, + {8, 11, 12}, + {10}, // remove + {13}, // remove + {1, 7}, + {3, 4, 6}, + {2, 7}, + {8, 11, 12}, + {0, 7}, + {4, 18, 20}, + {6, 16, 18}, + {3, 16, 20}, + {19}, // remove + {8}, // remove + {12}, // remove + {17, 21, 22, 23}, + {27}}; // remove +} + +BOOST_AUTO_TEST_CASE(zigzag_persistence_single) { + std::vector pairs; + auto stream = [&](int dim, ZP::Index birth, ZP::Index death) { pairs.emplace_back(dim, birth, death); }; + auto stream_inf = [&](int dim, ZP::Index birth) { pairs.emplace_back(dim, birth, -1); }; + ZP zp(stream, 28); + std::vector realIndices; + realIndices.reserve(13); + + std::vector > simplices = get_boundaries(); + + for (unsigned int i = 0; i < 14; ++i) { + zp.insert_cell(simplices[i], simplices[i].size() == 0 ? 0 : simplices[i].size() - 1); + } + + realIndices.emplace_back(0, 1, 3); + realIndices.emplace_back(0, 2, 4); + realIndices.emplace_back(0, 7, 8); + realIndices.emplace_back(1, 6, 10); + realIndices.emplace_back(0, 9, 11); + realIndices.emplace_back(1, 12, 13); + + for (unsigned int i = 14; i < 16; ++i) { + auto id = simplices[i][0]; + zp.remove_cell(id); + } + + for (unsigned int i = 16; i < 24; ++i) { + zp.insert_cell(simplices[i], simplices[i].size() == 0 ? 0 : simplices[i].size() - 1); + } + + realIndices.emplace_back(0, 5, 16); + realIndices.emplace_back(1, 14, 17); + realIndices.emplace_back(1, 15, 19); + realIndices.emplace_back(1, 20, 21); + realIndices.emplace_back(1, 18, 22); + + for (unsigned int i = 24; i < 27; ++i) { + auto id = simplices[i][0]; + zp.remove_cell(id); + } + + realIndices.emplace_back(1, 24, 25); + + zp.insert_cell(simplices[27], simplices[27].size() == 0 ? 0 : simplices[27].size() - 1); + + realIndices.emplace_back(2, 23, 27); + + auto id = simplices[28][0]; + zp.remove_cell(id); + + realIndices.emplace_back(0, 0, -1); + realIndices.emplace_back(0, 26, -1); + realIndices.emplace_back(2, 28, -1); + + auto start = pairs.size(); + zp.get_current_infinite_intervals(stream_inf); + std::sort(pairs.begin() + start, pairs.end(), [](const Interval& i1, const Interval& i2){ + if (i1.dim() != i2.dim()) return i1.dim() < i2.dim(); + return i1.birth() < i2.birth(); + }); + + test_indices(pairs, realIndices); +} + +BOOST_AUTO_TEST_CASE(zigzag_persistence_single_max1) { + std::vector pairs; + auto stream = [&](int dim, ZP::Index birth, ZP::Index death) { + if (dim < 1) pairs.emplace_back(dim, birth, death); + }; + auto stream_inf = [&](int dim, ZP::Index birth) { + if (dim < 1) pairs.emplace_back(dim, birth, -1); + }; + ZP zp(stream, 28); + std::vector realIndices; + realIndices.reserve(5); + + std::vector > simplices = get_boundaries(); + + for (unsigned int i = 0; i < 14; ++i) { + zp.insert_cell(simplices[i], simplices[i].size() == 0 ? 0 : simplices[i].size() - 1); + } + + realIndices.emplace_back(0, 1, 3); + realIndices.emplace_back(0, 2, 4); + realIndices.emplace_back(0, 7, 8); + realIndices.emplace_back(0, 9, 11); + + for (unsigned int i = 14; i < 16; ++i) { + auto id = simplices[i][0]; + zp.remove_cell(id); + } + + for (unsigned int i = 16; i < 24; ++i) { + zp.insert_cell(simplices[i], simplices[i].size() == 0 ? 0 : simplices[i].size() - 1); + } + + realIndices.emplace_back(0, 5, 16); + + for (unsigned int i = 24; i < 27; ++i) { + auto id = simplices[i][0]; + zp.remove_cell(id); + } + + zp.insert_cell(simplices[27], simplices[27].size() == 0 ? 0 : simplices[27].size() - 1); + auto id = simplices[28][0]; + zp.remove_cell(id); + + realIndices.emplace_back(0, 0, -1); + realIndices.emplace_back(0, 26, -1); + + auto start = pairs.size(); + zp.get_current_infinite_intervals(stream_inf); + std::sort(pairs.begin() + start, pairs.end(), [](const Interval& i1, const Interval& i2){ + if (i1.dim() != i2.dim()) return i1.dim() < i2.dim(); + return i1.birth() < i2.birth(); + }); + + test_indices(pairs, realIndices); +} diff --git a/src/common/doc/examples.h b/src/common/doc/examples.h index 2fa0c7c007..d363024beb 100644 --- a/src/common/doc/examples.h +++ b/src/common/doc/examples.h @@ -129,6 +129,12 @@ * @example example_one_skeleton_rips_from_distance_matrix.cpp * @example example_one_skeleton_rips_from_points.cpp * @example example_rips_complex_from_off_file.cpp + * \section Zigzag_persistence_example_section Zigzag_persistence + * @example example_usage_zigzag_persistence.cpp + * @example example_usage_filtered_zigzag_persistence.cpp + * @example example_usage_filtered_zigzag_persistence_with_storage.cpp + * @example example_zigzag_filtration_as_input_loop.cpp + * @example example_zzfiltration_from_file.cpp * \section Persistence_matrix_example_section Persistence_matrix * @example example_representative_cycles_from_matrix.cpp * @example example_simplex_tree_to_matrix.cpp diff --git a/src/common/doc/main_page.md b/src/common/doc/main_page.md index afd7734c41..334746b319 100644 --- a/src/common/doc/main_page.md +++ b/src/common/doc/main_page.md @@ -409,6 +409,31 @@ +### Zigzag Persistent Homology + + + + + + + + + + +
+ \image html "zigzag_ex.png" + + Zigzag filtrations are a generalization of (standard) filtrations: the inclusion maps can not only go forward, + but also backward. In other words, simplices can also be removed from the complex which is build. + The implemented algorithm is from \cite zigzag. + + Author: Clément Maria, Hannah Schreiber
+ Introduced in: GUDHI 3.11.0
+ Copyright: MIT
+
+ User manual: \ref zigzag_persistence +
+ ## Topological descriptors tools {#TopologicalDescriptorsTools} ### Bottleneck distance