diff --git a/flagser b/flagser index 7ffe927..da7f104 160000 --- a/flagser +++ b/flagser @@ -1 +1 @@ -Subproject commit 7ffe92782ae639b4265a52403d2d90d3151bb67b +Subproject commit da7f1048f444c5f9934a0a2bcab3dd31f9053395 diff --git a/pyflagser/flagser.py b/pyflagser/flagser.py index 934224f..f0c90b0 100644 --- a/pyflagser/flagser.py +++ b/pyflagser/flagser.py @@ -9,7 +9,8 @@ def flagser_unweighted(adjacency_matrix, min_dimension=0, max_dimension=np.inf, - directed=True, coeff=2, approximation=None): + directed=True, coeff=2, approximation=None, + in_memory=False): """Compute homology of a directed/undirected flag complex. From an adjacency_matrix construct all cells forming its associated flag @@ -48,6 +49,11 @@ def flagser_unweighted(adjacency_matrix, min_dimension=0, max_dimension=np.inf, computation. If ``None``, no approximation is made and all cells are used. For more details, please refer to [1]_. + in_memory: bool, optional, default: ``False`` + If ``True``, speeds up computation at the cost of increased memory + consumption, by enabling an efficient lookup in the data structure + instead of relying on hashing for looking up simplex indices. + Returns ------- out : dict of list @@ -101,7 +107,7 @@ def flagser_unweighted(adjacency_matrix, min_dimension=0, max_dimension=np.inf, # Call flagser binding homology = _compute_homology(vertices, edges, min_dimension, _max_dimension, directed, coeff, - _approximation, _filtration)[0] + _approximation, _filtration, in_memory)[0] # Creating dictionary of return values out = { @@ -114,7 +120,7 @@ def flagser_unweighted(adjacency_matrix, min_dimension=0, max_dimension=np.inf, def flagser_weighted(adjacency_matrix, max_edge_weight=None, min_dimension=0, max_dimension=np.inf, directed=True, filtration="max", - coeff=2, approximation=None): + coeff=2, approximation=None, in_memory=False): """Compute persistent homology of a directed/undirected filtered flag complex. @@ -184,6 +190,11 @@ def flagser_weighted(adjacency_matrix, max_edge_weight=None, min_dimension=0, computation. If ``None``, no approximation is made and all cells are used. For more details, please refer to [1]_. + in_memory: bool, optional, default: ``False`` + If ``True``, speeds up computation at the cost of increased memory + consumption, by enabling an efficient lookup in the data structure + instead of relying on hashing for looking up simplex indices. + Returns ------- out : dict of list @@ -247,7 +258,7 @@ def flagser_weighted(adjacency_matrix, max_edge_weight=None, min_dimension=0, # Call flagser binding homology = _compute_homology(vertices, edges, min_dimension, _max_dimension, directed, coeff, - _approximation, filtration)[0] + _approximation, filtration, in_memory)[0] # Create dictionary of return values out = { diff --git a/src/flagser_bindings.cpp b/src/flagser_bindings.cpp index 34cf1d6..eb3454a 100644 --- a/src/flagser_bindings.cpp +++ b/src/flagser_bindings.cpp @@ -1,11 +1,10 @@ -#include -#include +#include +#include #include #include - -#include -#include +#include +#include namespace py = pybind11; @@ -15,36 +14,79 @@ PYBIND11_MODULE(flagser_coeff_pybind, m) { PYBIND11_MODULE(flagser_pybind, m) { #endif + /* `persistence_computer_t` class requires to know if the computation is + * in memory or not. Because, `persistence_computer_t` contains the method + * to retrieve the different results of the homology computation. + * This class allows to store results from the computation independently if + * it was in memory or not. + * This duplicates the code in `persistence_computer_t`. + */ + class PersistenceValues { + public: + PersistenceValues(); + PersistenceValues(index_t eu, std::vector bn, + std::vector>> bd, + std::vector cc) + : euler_characteristic(eu), + betti_numbers(bn), + birth_deaths_by_dim(bd), + cell_count(cc) {} + index_t get_euler_characteristic() { return euler_characteristic; } + + std::vector get_betti_numbers() { return betti_numbers; } + + size_t get_betti_numbers(size_t dimension) { + return betti_numbers[dimension]; + } + + std::vector>> + get_persistence_diagram() { + return birth_deaths_by_dim; + } + + std::vector> get_persistence_diagram( + size_t dimension) { + return birth_deaths_by_dim[dimension]; + } + + std::vector get_cell_count() { return cell_count; } + + size_t get_cell_count(size_t dimension) { return cell_count[dimension]; } + + private: + index_t euler_characteristic = 0; + std::vector betti_numbers; + std::vector>> birth_deaths_by_dim; + std::vector cell_count; + }; + m.doc() = "Python interface for flagser"; m.attr("AVAILABLE_FILTRATIONS") = custom_filtration_computer; - using PersistenceComputer = - persistence_computer_t; - - py::class_(m, "PersistenceComputer", py::module_local()) + py::class_(m, "PersistenceValues", py::module_local()) .def("get_euler_characteristic", - &PersistenceComputer::get_euler_characteristic) + &PersistenceValues::get_euler_characteristic) .def("get_betti_numbers", - py::overload_cast<>(&PersistenceComputer::get_betti_numbers)) + py::overload_cast<>(&PersistenceValues::get_betti_numbers)) .def("get_betti_numbers", - py::overload_cast(&PersistenceComputer::get_betti_numbers)) + py::overload_cast(&PersistenceValues::get_betti_numbers)) .def("get_cell_count", - py::overload_cast<>(&PersistenceComputer::get_cell_count)) + py::overload_cast<>(&PersistenceValues::get_cell_count)) .def("get_cell_count", - py::overload_cast(&PersistenceComputer::get_cell_count)) + py::overload_cast(&PersistenceValues::get_cell_count)) .def("get_persistence_diagram", - py::overload_cast<>(&PersistenceComputer::get_persistence_diagram)) + py::overload_cast<>(&PersistenceValues::get_persistence_diagram)) .def("get_persistence_diagram", py::overload_cast( - &PersistenceComputer::get_persistence_diagram)); + &PersistenceValues::get_persistence_diagram)); m.def("compute_homology", [](std::vector& vertices, std::vector>& edges, unsigned short min_dim, short max_dim, bool directed, coefficient_t modulus, - signed int approximation, - std::string filtration) { + signed int approximation, std::string filtration, + bool in_memory) { flagser_parameters params; // Save std::cout status @@ -79,6 +121,10 @@ PYBIND11_MODULE(flagser_pybind, m) { // Calls Trivial output, disable the generation of an output file params.output_format = std::string("none"); + // Enable in memory computation for directed flags + // Threading memory for performances + params.input_format = in_memory; + // Building the filtered directed graph auto graph = filtered_directed_graph_t(vertices, params.directed); @@ -117,12 +163,37 @@ PYBIND11_MODULE(flagser_pybind, m) { // Disable cout std::cout.rdbuf(nullptr); + // Will contain Homology results + std::vector vec_pv; + // Running flagser's compute_homology routine - auto subgraph_persistence_computer = compute_homology(graph, params); + if (params.in_memory) { + auto subgraph_persistence_computer = + compute_homology( + graph, params); + /* FIXME: right now, the loop needs to be duplicated because in order + * right now, it's not possible to create a vector of optional's + */ + for (auto& sub : subgraph_persistence_computer) + vec_pv.emplace_back(PersistenceValues( + sub.get_euler_characteristic(), sub.get_betti_numbers(), + sub.get_persistence_diagram(), sub.get_cell_count())); + + } else { + auto subgraph_persistence_computer = compute_homology< + directed_flag_complex_computer::directed_flag_complex_computer_t>( + graph, params); + /* FIXME: Remove when supported the duplicated loop */ + for (auto& sub : subgraph_persistence_computer) + vec_pv.emplace_back(PersistenceValues( + sub.get_euler_characteristic(), sub.get_betti_numbers(), + sub.get_persistence_diagram(), sub.get_cell_count())); + } // Re-enable again cout std::cout.rdbuf(cout_buff); - return subgraph_persistence_computer; + return vec_pv; }); }