diff --git a/src/Persistence_matrix/include/gudhi/Persistence_matrix/Chain_matrix.h b/src/Persistence_matrix/include/gudhi/Persistence_matrix/Chain_matrix.h index 7611990b8..dce0a1d50 100644 --- a/src/Persistence_matrix/include/gudhi/Persistence_matrix/Chain_matrix.h +++ b/src/Persistence_matrix/include/gudhi/Persistence_matrix/Chain_matrix.h @@ -913,24 +913,28 @@ inline void Chain_matrix::print() const { std::cout << "Column Matrix:\n"; if constexpr (!Master_matrix::Option_list::has_map_column_container) { - for (ID_index i = 0; i < pivotToColumnIndex_.size() && pivotToColumnIndex_[i] != static_cast(-1); ++i) { + for (ID_index i = 0; i < pivotToColumnIndex_.size(); ++i) { Index pos = pivotToColumnIndex_[i]; - const Column& col = matrix_[pos]; - for (const auto& entry : col) { - std::cout << entry.get_row_index() << " "; + if (pos != static_cast(-1)){ + const Column& col = matrix_[pos]; + for (const auto& entry : col) { + std::cout << entry.get_row_index() << " "; + } + std::cout << "(" << i << ", " << pos << ")\n"; } - std::cout << "(" << i << ", " << pos << ")\n"; } if constexpr (Master_matrix::Option_list::has_row_access) { std::cout << "\n"; std::cout << "Row Matrix:\n"; - for (ID_index i = 0; i < pivotToColumnIndex_.size() && pivotToColumnIndex_[i] != static_cast(-1); ++i) { + for (ID_index i = 0; i < pivotToColumnIndex_.size(); ++i) { Index pos = pivotToColumnIndex_[i]; - const Row& row = RA_opt::get_row(pos); - for (const auto& entry : row) { - std::cout << entry.get_column_index() << " "; + if (pos != static_cast(-1)){ + const Row& row = RA_opt::get_row(pos); + for (const auto& entry : row) { + std::cout << entry.get_column_index() << " "; + } + std::cout << "(" << i << ", " << pos << ")\n"; } - std::cout << "(" << i << ", " << pos << ")\n"; } } } else { diff --git a/src/Persistence_matrix/include/gudhi/Persistence_matrix/base_pairing.h b/src/Persistence_matrix/include/gudhi/Persistence_matrix/base_pairing.h index 8ce95ce74..780040f01 100644 --- a/src/Persistence_matrix/include/gudhi/Persistence_matrix/base_pairing.h +++ b/src/Persistence_matrix/include/gudhi/Persistence_matrix/base_pairing.h @@ -132,10 +132,11 @@ inline const typename Base_pairing::Barcode& Base_pairing inline void Base_pairing::_reduce() { - std::unordered_map pivotsToColumn(_matrix()->get_number_of_columns()); + std::unordered_map negativeColumns(_matrix()->get_number_of_columns()); auto dim = _matrix()->get_max_dimension(); std::vector > columnsByDim(dim + 1); + for (auto& v : columnsByDim) v.reserve(_matrix()->get_number_of_columns()); for (unsigned int i = 0; i < _matrix()->get_number_of_columns(); i++) { columnsByDim[dim - _matrix()->get_column_dimension(i)].push_back(i); } @@ -144,17 +145,21 @@ inline void Base_pairing::_reduce() for (Index i : cols) { auto& curr = _matrix()->get_column(i); if (curr.is_empty()) { - if (pivotsToColumn.find(i) == pivotsToColumn.end()) { + if (negativeColumns.find(i) == negativeColumns.end()) { barcode_.emplace_back(i, -1, dim); } } else { ID_index pivot = curr.get_pivot(); + auto it = idToPosition_.find(pivot); + Index pivotColumnNumber = it == idToPosition_.end() ? pivot : it->second; + auto itNeg = negativeColumns.find(pivotColumnNumber); + Index pivotKiller = itNeg == negativeColumns.end() ? -1 : itNeg->second; - while (pivot != static_cast(-1) && pivotsToColumn.find(pivot) != pivotsToColumn.end()) { + while (pivot != static_cast(-1) && pivotKiller != static_cast(-1)) { if constexpr (Master_matrix::Option_list::is_z2) { - curr += _matrix()->get_column(pivotsToColumn.at(pivot)); + curr += _matrix()->get_column(pivotKiller); } else { - auto& toadd = _matrix()->get_column(pivotsToColumn.at(pivot)); + auto& toadd = _matrix()->get_column(pivotKiller); typename Master_matrix::Element coef = toadd.get_pivot_value(); auto& operators = _matrix()->colSettings_->operators; coef = operators.get_inverse(coef); @@ -163,12 +168,14 @@ inline void Base_pairing::_reduce() } pivot = curr.get_pivot(); + it = idToPosition_.find(pivot); + pivotColumnNumber = it == idToPosition_.end() ? pivot : it->second; + itNeg = negativeColumns.find(pivotColumnNumber); + pivotKiller = itNeg == negativeColumns.end() ? -1 : itNeg->second; } if (pivot != static_cast(-1)) { - pivotsToColumn.emplace(pivot, i); - auto it = idToPosition_.find(pivot); - auto pivotColumnNumber = it == idToPosition_.end() ? pivot : it->second; + negativeColumns.emplace(pivotColumnNumber, i); _matrix()->get_column(pivotColumnNumber).clear(); barcode_.emplace_back(pivotColumnNumber, i, dim - 1); } else { diff --git a/src/Persistence_matrix/test/pm_matrix_tests.h b/src/Persistence_matrix/test/pm_matrix_tests.h index 3d0aac966..0028079c5 100644 --- a/src/Persistence_matrix/test/pm_matrix_tests.h +++ b/src/Persistence_matrix/test/pm_matrix_tests.h @@ -1348,8 +1348,9 @@ template void test_barcode() { struct BarComp { bool operator()(const std::tuple& c1, const std::tuple& c2) const { - if (std::get<0>(c1) == std::get<0>(c2)) return std::get<1>(c1) < std::get<1>(c2); - return std::get<0>(c1) < std::get<0>(c2); + if (std::get<0>(c1) != std::get<0>(c2)) return std::get<0>(c1) < std::get<0>(c2); + if (std::get<1>(c1) != std::get<1>(c2)) return std::get<1>(c1) < std::get<1>(c2); + return std::get<2>(c1) < std::get<2>(c2); } }; @@ -1420,12 +1421,220 @@ void test_barcode() { } template -void test_shifted_barcode() { +void test_shifted_barcode1() { + using C = typename Matrix::Column; + struct BarComp { + bool operator()(const std::tuple& c1, const std::tuple& c2) const { + if (std::get<0>(c1) != std::get<0>(c2)) return std::get<0>(c1) < std::get<0>(c2); + if (std::get<1>(c1) != std::get<1>(c2)) return std::get<1>(c1) < std::get<1>(c2); + return std::get<2>(c1) < std::get<2>(c2); + } + }; + + Matrix m(17, 2); + if constexpr (is_z2()) { + m.insert_boundary(0, {}, 0); + m.insert_boundary(1, {}, 0); + m.insert_boundary(2, {}, 0); + m.insert_boundary(3, {}, 0); + m.insert_boundary(4, {}, 0); + m.insert_boundary(5, {}, 0); + m.insert_boundary(6, {}, 0); + m.insert_boundary(10, {0, 1}, 1); + m.insert_boundary(11, {1, 3}, 1); + m.insert_boundary(12, {2, 3}, 1); + m.insert_boundary(13, {2, 4}, 1); + m.insert_boundary(14, {3, 4}, 1); + m.insert_boundary(15, {2, 6}, 1); + m.insert_boundary(16, {4, 6}, 1); + m.insert_boundary(17, {5, 6}, 1); + m.insert_boundary(30, {12, 13, 14}, 2); + m.insert_boundary(31, {13, 15, 16}, 2); + } else { + m.insert_boundary(0, {}, 0); + m.insert_boundary(1, {}, 0); + m.insert_boundary(2, {}, 0); + m.insert_boundary(3, {}, 0); + m.insert_boundary(4, {}, 0); + m.insert_boundary(5, {}, 0); + m.insert_boundary(6, {}, 0); + m.insert_boundary(10, {{0, 1}, {1, 1}}, 1); + m.insert_boundary(11, {{1, 1}, {3, 1}}, 1); + m.insert_boundary(12, {{2, 1}, {3, 1}}, 1); + m.insert_boundary(13, {{2, 1}, {4, 1}}, 1); + m.insert_boundary(14, {{3, 1}, {4, 1}}, 1); + m.insert_boundary(15, {{2, 1}, {6, 1}}, 1); + m.insert_boundary(16, {{4, 1}, {6, 1}}, 1); + m.insert_boundary(17, {{5, 1}, {6, 1}}, 1); + m.insert_boundary(30, {{12, 1}, {13, 1}, {14, 1}}, 2); + m.insert_boundary(31, {{13, 1}, {15, 1}, {16, 1}}, 2); + } + + const auto& barcode = m.get_current_barcode(); + + std::vector > reducedMatrix; + if constexpr (is_z2()) { + if constexpr (Matrix::Option_list::is_of_boundary_type) { + reducedMatrix.emplace_back(); + reducedMatrix.emplace_back(); + reducedMatrix.emplace_back(); + reducedMatrix.emplace_back(); + reducedMatrix.emplace_back(); + reducedMatrix.emplace_back(); + reducedMatrix.emplace_back(); + reducedMatrix.push_back({0, 1}); + reducedMatrix.push_back({1, 3}); + reducedMatrix.push_back({1, 2}); + reducedMatrix.push_back({2, 4}); + reducedMatrix.emplace_back(); + reducedMatrix.push_back({2, 6}); + reducedMatrix.emplace_back(); + reducedMatrix.push_back({2, 5}); + reducedMatrix.push_back({12, 13, 14}); + reducedMatrix.push_back({13, 15, 16}); + } else { + reducedMatrix.push_back({0}); + reducedMatrix.push_back({0, 1}); + reducedMatrix.push_back({0, 2}); + reducedMatrix.push_back({0, 3}); + reducedMatrix.push_back({0, 4}); + reducedMatrix.push_back({0, 5}); + reducedMatrix.push_back({0, 6}); + reducedMatrix.push_back({10}); + reducedMatrix.push_back({10, 11}); + reducedMatrix.push_back({10, 11, 12}); + reducedMatrix.push_back({10, 11, 12, 13}); + reducedMatrix.push_back({12, 13, 14}); + reducedMatrix.push_back({10, 11, 12, 15}); + reducedMatrix.push_back({13, 15, 16}); + reducedMatrix.push_back({10, 11, 12, 15, 17}); + reducedMatrix.push_back({30}); + reducedMatrix.push_back({31}); + } + } else { + if constexpr (Matrix::Option_list::is_of_boundary_type) { + reducedMatrix.emplace_back(); + reducedMatrix.emplace_back(); + reducedMatrix.emplace_back(); + reducedMatrix.emplace_back(); + reducedMatrix.emplace_back(); + reducedMatrix.emplace_back(); + reducedMatrix.emplace_back(); + reducedMatrix.push_back({{0, 1}, {1, 1}}); + reducedMatrix.push_back({{1, 1}, {3, 1}}); + reducedMatrix.push_back({{1, 1}, {2, 1}}); + reducedMatrix.push_back({{2, 1}, {4, 1}}); + reducedMatrix.emplace_back(); + reducedMatrix.push_back({{2, 1}, {6, 1}}); + reducedMatrix.emplace_back(); + reducedMatrix.push_back({{2, 1}, {5, 1}}); + reducedMatrix.push_back({{12, 1}, {13, 1}, {14, 1}}); + reducedMatrix.push_back({{13, 1}, {15, 1}, {16, 1}}); + } else { + reducedMatrix.push_back({{0, 1}}); + reducedMatrix.push_back({{0, 1}, {1, 1}}); + reducedMatrix.push_back({{0, 1}, {2, 1}}); + reducedMatrix.push_back({{0, 1}, {3, 1}}); + reducedMatrix.push_back({{0, 1}, {4, 1}}); + reducedMatrix.push_back({{0, 1}, {5, 1}}); + reducedMatrix.push_back({{0, 1}, {6, 1}}); + reducedMatrix.push_back({{10, 1}}); + reducedMatrix.push_back({{10, 1}, {11, 1}}); + reducedMatrix.push_back({{10, 1}, {11, 1}, {12, 1}}); + reducedMatrix.push_back({{10, 1}, {11, 1}, {12, 1}, {13, 1}}); + reducedMatrix.push_back({{12, 1}, {13, 1}, {14, 1}}); + reducedMatrix.push_back({{10, 1}, {11, 1}, {12, 1}, {15, 1}}); + reducedMatrix.push_back({{13, 1}, {15, 1}, {16, 1}}); + reducedMatrix.push_back({{10, 1}, {11, 1}, {12, 1}, {15, 1}, {17, 1}}); + reducedMatrix.push_back({{30, 1}}); + reducedMatrix.push_back({{31, 1}}); + } + } + + if constexpr (Matrix::Option_list::column_indexation_type == Column_indexation_types::IDENTIFIER){ + test_column_equality(reducedMatrix[0], get_column_content_via_iterators(m.get_column(0))); + test_column_equality(reducedMatrix[1], get_column_content_via_iterators(m.get_column(1))); + test_column_equality(reducedMatrix[2], get_column_content_via_iterators(m.get_column(2))); + test_column_equality(reducedMatrix[3], get_column_content_via_iterators(m.get_column(3))); + test_column_equality(reducedMatrix[4], get_column_content_via_iterators(m.get_column(4))); + test_column_equality(reducedMatrix[5], get_column_content_via_iterators(m.get_column(5))); + test_column_equality(reducedMatrix[6], get_column_content_via_iterators(m.get_column(6))); + test_column_equality(reducedMatrix[7], get_column_content_via_iterators(m.get_column(10))); + test_column_equality(reducedMatrix[8], get_column_content_via_iterators(m.get_column(11))); + test_column_equality(reducedMatrix[9], get_column_content_via_iterators(m.get_column(12))); + test_column_equality(reducedMatrix[10], get_column_content_via_iterators(m.get_column(13))); + test_column_equality(reducedMatrix[11], get_column_content_via_iterators(m.get_column(14))); + test_column_equality(reducedMatrix[12], get_column_content_via_iterators(m.get_column(15))); + test_column_equality(reducedMatrix[13], get_column_content_via_iterators(m.get_column(16))); + test_column_equality(reducedMatrix[14], get_column_content_via_iterators(m.get_column(17))); + test_column_equality(reducedMatrix[15], get_column_content_via_iterators(m.get_column(30))); + test_column_equality(reducedMatrix[16], get_column_content_via_iterators(m.get_column(31))); + } else { + test_content_equality(reducedMatrix, m); + } + + std::set, BarComp> bars1; + std::set, BarComp> bars2; + std::set, BarComp> bars3; + // bars are not ordered the same for all matrices + for (auto it = barcode.begin(); it != barcode.end(); ++it) { + //three access possibilities + bars1.emplace(it->dim, it->birth, it->death); + bars2.emplace(std::get<2>(*it), std::get<0>(*it), std::get<1>(*it)); + auto [ x, y, z ] = *it; + bars3.emplace(z, x, y); + } + auto it = bars1.begin(); + BOOST_CHECK_EQUAL(std::get<0>(*it), 0); + BOOST_CHECK_EQUAL(std::get<1>(*it), 0); + BOOST_CHECK_EQUAL(std::get<2>(*it), -1); + ++it; + BOOST_CHECK_EQUAL(std::get<0>(*it), 0); + BOOST_CHECK_EQUAL(std::get<1>(*it), 1); + BOOST_CHECK_EQUAL(std::get<2>(*it), 7); + ++it; + BOOST_CHECK_EQUAL(std::get<0>(*it), 0); + BOOST_CHECK_EQUAL(std::get<1>(*it), 2); + BOOST_CHECK_EQUAL(std::get<2>(*it), 9); + ++it; + BOOST_CHECK_EQUAL(std::get<0>(*it), 0); + BOOST_CHECK_EQUAL(std::get<1>(*it), 3); + BOOST_CHECK_EQUAL(std::get<2>(*it), 8); + ++it; + BOOST_CHECK_EQUAL(std::get<0>(*it), 0); + BOOST_CHECK_EQUAL(std::get<1>(*it), 4); + BOOST_CHECK_EQUAL(std::get<2>(*it), 10); + ++it; + BOOST_CHECK_EQUAL(std::get<0>(*it), 0); + BOOST_CHECK_EQUAL(std::get<1>(*it), 5); + BOOST_CHECK_EQUAL(std::get<2>(*it), 14); + ++it; + BOOST_CHECK_EQUAL(std::get<0>(*it), 0); + BOOST_CHECK_EQUAL(std::get<1>(*it), 6); + BOOST_CHECK_EQUAL(std::get<2>(*it), 12); + ++it; + BOOST_CHECK_EQUAL(std::get<0>(*it), 1); + BOOST_CHECK_EQUAL(std::get<1>(*it), 11); + BOOST_CHECK_EQUAL(std::get<2>(*it), 15); + ++it; + BOOST_CHECK_EQUAL(std::get<0>(*it), 1); + BOOST_CHECK_EQUAL(std::get<1>(*it), 13); + BOOST_CHECK_EQUAL(std::get<2>(*it), 16); + ++it; + BOOST_CHECK(it == bars1.end()); + + BOOST_CHECK(bars1 == bars2); + BOOST_CHECK(bars1 == bars3); +} + +template +void test_shifted_barcode2() { using C = typename Matrix::Column; struct BarComp { bool operator()(const std::tuple& c1, const std::tuple& c2) const { - if (std::get<0>(c1) == std::get<0>(c2)) return std::get<1>(c1) < std::get<1>(c2); - return std::get<0>(c1) < std::get<0>(c2); + if (std::get<0>(c1) != std::get<0>(c2)) return std::get<0>(c1) < std::get<0>(c2); + if (std::get<1>(c1) != std::get<1>(c2)) return std::get<1>(c1) < std::get<1>(c2); + return std::get<2>(c1) < std::get<2>(c2); } }; @@ -1553,6 +1762,12 @@ void test_shifted_barcode() { BOOST_CHECK(bars1 == bars3); } +template +void test_shifted_barcode() { + test_shifted_barcode1(); + test_shifted_barcode2(); +} + template void test_base_swaps() { auto columns = build_simple_boundary_matrix();