All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
Chain_matrix.h
Go to the documentation of this file.
1/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
2 * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
3 * Author(s): Hannah Schreiber
4 *
5 * Copyright (C) 2022-24 Inria
6 *
7 * Modification(s):
8 * - YYYY/MM Author: Description of the modification
9 */
10
17#ifndef PM_CHAIN_MATRIX_H
18#define PM_CHAIN_MATRIX_H
19
20#include <iostream> //print() only
21#include <set>
22#include <map>
23#include <stdexcept>
24#include <vector>
25#include <utility> //std::swap, std::move & std::exchange
26#include <algorithm> //std::sort
27
29
30namespace Gudhi {
31namespace persistence_matrix {
32
44template <class Master_matrix>
45class Chain_matrix : public Master_matrix::Matrix_dimension_option,
46 public Master_matrix::Chain_pairing_option,
47 public Master_matrix::Chain_vine_swap_option,
48 public Master_matrix::Chain_representative_cycles_option,
49 public Master_matrix::Matrix_row_access_option
50{
51 public:
55 using Field_operators = typename Master_matrix::Field_operators;
56 using Field_element = typename Master_matrix::Element;
57 using Column = typename Master_matrix::Column;
58 using Row = typename Master_matrix::Row;
60 using Entry = typename Master_matrix::Matrix_entry;
61 using Entry_constructor = typename Master_matrix::Entry_constructor;
62 using Column_settings = typename Master_matrix::Column_settings;
64 using Boundary = typename Master_matrix::Boundary;
65 using Entry_representative = typename Master_matrix::Entry_representative;
66 using Index = typename Master_matrix::Index;
67 using ID_index = typename Master_matrix::ID_index;
68 using Pos_index = typename Master_matrix::Pos_index;
69 using Dimension = typename Master_matrix::Dimension;
79 Chain_matrix(Column_settings* colSettings);
104 template <class Boundary_range = Boundary>
105 Chain_matrix(const std::vector<Boundary_range>& orderedBoundaries,
106 Column_settings* colSettings);
116 Chain_matrix(unsigned int numberOfColumns, Column_settings* colSettings);
138 template <typename BirthComparatorFunction, typename DeathComparatorFunction>
139 Chain_matrix(Column_settings* colSettings,
140 const BirthComparatorFunction& birthComparator,
141 const DeathComparatorFunction& deathComparator);
178 template <typename BirthComparatorFunction, typename DeathComparatorFunction, class Boundary_range = Boundary>
179 Chain_matrix(const std::vector<Boundary_range>& orderedBoundaries,
180 Column_settings* colSettings,
181 const BirthComparatorFunction& birthComparator,
182 const DeathComparatorFunction& deathComparator);
205 template <typename BirthComparatorFunction, typename DeathComparatorFunction>
206 Chain_matrix(unsigned int numberOfColumns,
207 Column_settings* colSettings,
208 const BirthComparatorFunction& birthComparator,
209 const DeathComparatorFunction& deathComparator);
219 Chain_matrix(const Chain_matrix& matrixToCopy,
220 Column_settings* colSettings = nullptr);
226 Chain_matrix(Chain_matrix&& other) noexcept;
227
249 template <class Boundary_range = Boundary>
250 std::vector<Entry_representative> insert_boundary(
251 const Boundary_range& boundary, Dimension dim = Master_matrix::template get_null_value<Dimension>());
270 template <class Boundary_range = Boundary>
271 std::vector<Entry_representative> insert_boundary(
272 ID_index cellID,
273 const Boundary_range& boundary,
274 Dimension dim = Master_matrix::template get_null_value<Dimension>());
282 Column& get_column(Index columnIndex);
290 const Column& get_column(Index columnIndex) const;
307 void remove_maximal_cell(ID_index cellID);
331 void remove_maximal_cell(ID_index cellID, const std::vector<ID_index>& columnsToSwap);
346 void remove_last();
347
354
361 Dimension get_column_dimension(Index columnIndex) const;
362
373 void add_to(Index sourceColumnIndex, Index targetColumnIndex);
386 void multiply_target_and_add_to(Index sourceColumnIndex,
387 const Field_element& coefficient,
388 Index targetColumnIndex);
401 void multiply_source_and_add_to(const Field_element& coefficient,
402 Index sourceColumnIndex,
403 Index targetColumnIndex);
404
413 bool is_zero_entry(Index columnIndex, ID_index rowIndex) const;
422 bool is_zero_column(Index columnIndex);
423
437 ID_index get_pivot(Index columnIndex);
438
445 void reset(Column_settings* colSettings) {
446 matrix_.clear();
447 pivotToColumnIndex_.clear();
448 nextIndex_ = 0;
449 colSettings_ = colSettings;
450 }
451
455 Chain_matrix& operator=(const Chain_matrix& other);
459 friend void swap(Chain_matrix& matrix1, Chain_matrix& matrix2) {
460 swap(static_cast<typename Master_matrix::Matrix_dimension_option&>(matrix1),
461 static_cast<typename Master_matrix::Matrix_dimension_option&>(matrix2));
462 swap(static_cast<typename Master_matrix::Chain_pairing_option&>(matrix1),
463 static_cast<typename Master_matrix::Chain_pairing_option&>(matrix2));
464 swap(static_cast<typename Master_matrix::Chain_vine_swap_option&>(matrix1),
465 static_cast<typename Master_matrix::Chain_vine_swap_option&>(matrix2));
466 swap(static_cast<typename Master_matrix::Chain_representative_cycles_option&>(matrix1),
467 static_cast<typename Master_matrix::Chain_representative_cycles_option&>(matrix2));
468 matrix1.matrix_.swap(matrix2.matrix_);
469 matrix1.pivotToColumnIndex_.swap(matrix2.pivotToColumnIndex_);
470 std::swap(matrix1.nextIndex_, matrix2.nextIndex_);
471 std::swap(matrix1.colSettings_, matrix2.colSettings_);
472
474 swap(static_cast<typename Master_matrix::Matrix_row_access_option&>(matrix1),
475 static_cast<typename Master_matrix::Matrix_row_access_option&>(matrix2));
476 }
477 }
478
479 void print() const; // for debug
480
481 friend class Id_to_index_overlay<Chain_matrix<Master_matrix>, Master_matrix>;
482
483 private:
484 using Dim_opt = typename Master_matrix::Matrix_dimension_option;
485 using Swap_opt = typename Master_matrix::Chain_vine_swap_option;
486 using Pair_opt = typename Master_matrix::Chain_pairing_option;
487 using Rep_opt = typename Master_matrix::Chain_representative_cycles_option;
488 using RA_opt = typename Master_matrix::Matrix_row_access_option;
489 using Column_container = typename Master_matrix::Column_container;
490 using Dictionary = typename Master_matrix::template Dictionary<Index>;
491 using Barcode = typename Master_matrix::Barcode;
492 using Bar_dictionary = typename Master_matrix::Bar_dictionary;
493 using Tmp_column = typename std::conditional<Master_matrix::Option_list::is_z2,
494 std::set<ID_index>,
495 std::map<ID_index, Field_element>
496 >::type;
497
498 Column_container matrix_;
499 Dictionary pivotToColumnIndex_;
500 Index nextIndex_;
501 Column_settings* colSettings_;
503 template <class Boundary_range>
504 std::vector<Entry_representative> _reduce_boundary(ID_index cellID, const Boundary_range& boundary, Dimension dim);
505 void _reduce_by_G(Tmp_column& column, std::vector<Entry_representative>& chainsInH, Index currentPivot);
506 void _reduce_by_F(Tmp_column& column, std::vector<Entry_representative>& chainsInF, Index currentPivot);
507 void _build_from_H(ID_index cellID, Tmp_column& column, std::vector<Entry_representative>& chainsInH);
508 void _update_largest_death_in_F(const std::vector<Entry_representative>& chainsInF);
509 void _insert_chain(const Tmp_column& column, Dimension dimension);
510 void _insert_chain(const Tmp_column& column, Dimension dimension, Index pair);
511 void _add_to(const Column& column, Tmp_column& set, unsigned int coef);
512 template <typename F>
513 void _add_to(Column& target, F&& addition);
514 void _remove_last(Index lastIndex);
515 void _update_barcode(Pos_index birth);
516 void _add_bar(Dimension dim);
517 template <class Container>
518 void _container_insert(const Container& column, Index pos, Dimension dim);
519 void _container_insert(const Column& column, [[maybe_unused]] Index pos = 0);
520
521 constexpr Barcode& _barcode();
522 constexpr Bar_dictionary& _indexToBar();
523 constexpr Pos_index& _nextPosition();
524};
525
526template <class Master_matrix>
528 : Dim_opt(Master_matrix::template get_null_value<Dimension>()),
529 Pair_opt(),
530 Swap_opt(),
531 Rep_opt(),
532 RA_opt(),
533 nextIndex_(0),
534 colSettings_(colSettings)
535{}
536
537template <class Master_matrix>
538template <class Boundary_range>
539inline Chain_matrix<Master_matrix>::Chain_matrix(const std::vector<Boundary_range>& orderedBoundaries,
540 Column_settings* colSettings)
541 : Dim_opt(Master_matrix::template get_null_value<Dimension>()),
542 Pair_opt(),
543 Swap_opt(),
544 Rep_opt(),
545 RA_opt(orderedBoundaries.size()),
546 nextIndex_(0),
547 colSettings_(colSettings)
548{
549 matrix_.reserve(orderedBoundaries.size());
551 pivotToColumnIndex_.reserve(orderedBoundaries.size());
552 } else {
553 pivotToColumnIndex_.resize(orderedBoundaries.size(), Master_matrix::template get_null_value<Index>());
554 }
555
556 for (const Boundary_range& b : orderedBoundaries) {
558 }
559}
560
561template <class Master_matrix>
562inline Chain_matrix<Master_matrix>::Chain_matrix(unsigned int numberOfColumns,
563 Column_settings* colSettings)
564 : Dim_opt(Master_matrix::template get_null_value<Dimension>()),
565 Pair_opt(),
566 Swap_opt(),
567 Rep_opt(),
568 RA_opt(numberOfColumns),
569 nextIndex_(0),
570 colSettings_(colSettings)
571{
572 matrix_.reserve(numberOfColumns);
574 pivotToColumnIndex_.reserve(numberOfColumns);
575 } else {
576 pivotToColumnIndex_.resize(numberOfColumns, Master_matrix::template get_null_value<Index>());
577 }
578}
579
580template <class Master_matrix>
581template <typename BirthComparatorFunction, typename DeathComparatorFunction>
583 const BirthComparatorFunction& birthComparator,
584 const DeathComparatorFunction& deathComparator)
585 : Dim_opt(Master_matrix::template get_null_value<Dimension>()),
586 Pair_opt(),
587 Swap_opt(birthComparator, deathComparator),
588 Rep_opt(),
589 RA_opt(),
590 nextIndex_(0),
591 colSettings_(colSettings)
592{}
593
594template <class Master_matrix>
595template <typename BirthComparatorFunction, typename DeathComparatorFunction, class Boundary_range>
596inline Chain_matrix<Master_matrix>::Chain_matrix(const std::vector<Boundary_range>& orderedBoundaries,
597 Column_settings* colSettings,
598 const BirthComparatorFunction& birthComparator,
599 const DeathComparatorFunction& deathComparator)
600 : Dim_opt(Master_matrix::template get_null_value<Dimension>()),
601 Pair_opt(),
602 Swap_opt(birthComparator, deathComparator),
603 Rep_opt(),
604 RA_opt(orderedBoundaries.size()),
605 nextIndex_(0),
606 colSettings_(colSettings)
607{
608 matrix_.reserve(orderedBoundaries.size());
610 pivotToColumnIndex_.reserve(orderedBoundaries.size());
611 } else {
612 pivotToColumnIndex_.resize(orderedBoundaries.size(), Master_matrix::template get_null_value<Index>());
613 }
614 for (const Boundary_range& b : orderedBoundaries) {
616 }
617}
618
619template <class Master_matrix>
620template <typename BirthComparatorFunction, typename DeathComparatorFunction>
621inline Chain_matrix<Master_matrix>::Chain_matrix(unsigned int numberOfColumns,
622 Column_settings* colSettings,
623 const BirthComparatorFunction& birthComparator,
624 const DeathComparatorFunction& deathComparator)
625 : Dim_opt(Master_matrix::template get_null_value<Dimension>()),
626 Pair_opt(),
627 Swap_opt(birthComparator, deathComparator),
628 Rep_opt(),
629 RA_opt(numberOfColumns),
630 nextIndex_(0),
631 colSettings_(colSettings)
632{
633 matrix_.reserve(numberOfColumns);
635 pivotToColumnIndex_.reserve(numberOfColumns);
636 } else {
637 pivotToColumnIndex_.resize(numberOfColumns, Master_matrix::template get_null_value<Index>());
638 }
639}
640
641template <class Master_matrix>
643 : Dim_opt(static_cast<const Dim_opt&>(matrixToCopy)),
644 Pair_opt(static_cast<const Pair_opt&>(matrixToCopy)),
645 Swap_opt(static_cast<const Swap_opt&>(matrixToCopy)),
646 Rep_opt(static_cast<const Rep_opt&>(matrixToCopy)),
647 RA_opt(static_cast<const RA_opt&>(matrixToCopy)),
648 pivotToColumnIndex_(matrixToCopy.pivotToColumnIndex_),
649 nextIndex_(matrixToCopy.nextIndex_),
650 colSettings_(colSettings == nullptr ? matrixToCopy.colSettings_ : colSettings)
651{
652 matrix_.reserve(matrixToCopy.matrix_.size());
653 for (const auto& cont : matrixToCopy.matrix_){
655 _container_insert(cont.second, cont.first);
656 } else {
657 _container_insert(cont);
658 }
659 }
660}
661
662template <class Master_matrix>
664 : Dim_opt(std::move(static_cast<Dim_opt&>(other))),
665 Pair_opt(std::move(static_cast<Pair_opt&>(other))),
666 Swap_opt(std::move(static_cast<Swap_opt&>(other))),
667 Rep_opt(std::move(static_cast<Rep_opt&>(other))),
668 RA_opt(std::move(static_cast<RA_opt&>(other))),
669 matrix_(std::move(other.matrix_)),
670 pivotToColumnIndex_(std::move(other.pivotToColumnIndex_)),
671 nextIndex_(std::exchange(other.nextIndex_, 0)),
672 colSettings_(std::exchange(other.colSettings_, nullptr))
673{}
674
675template <class Master_matrix>
676template <class Boundary_range>
677inline std::vector<typename Master_matrix::Entry_representative> Chain_matrix<Master_matrix>::insert_boundary(
678 const Boundary_range& boundary, Dimension dim)
679{
680 return insert_boundary(nextIndex_, boundary, dim);
681}
682
683template <class Master_matrix>
684template <class Boundary_range>
685inline std::vector<typename Master_matrix::Entry_representative> Chain_matrix<Master_matrix>::insert_boundary(
686 ID_index cellID, const Boundary_range& boundary, Dimension dim)
687{
688 if constexpr (!Master_matrix::Option_list::has_map_column_container) {
689 if (pivotToColumnIndex_.size() <= cellID) {
690 pivotToColumnIndex_.resize(cellID * 2 + 1, Master_matrix::template get_null_value<Index>());
691 }
692 }
693
694 if constexpr (Master_matrix::Option_list::has_vine_update && Master_matrix::Option_list::has_column_pairings) {
695 if constexpr (Master_matrix::Option_list::has_map_column_container) {
696 Swap_opt::CP::pivotToPosition_.try_emplace(cellID, _nextPosition());
697 } else {
698 if (Swap_opt::CP::pivotToPosition_.size() <= cellID)
699 Swap_opt::CP::pivotToPosition_.resize(pivotToColumnIndex_.size(),
700 Master_matrix::template get_null_value<Pos_index>());
701 Swap_opt::CP::pivotToPosition_[cellID] = _nextPosition();
702 }
703 }
704
705 if constexpr (Master_matrix::Option_list::has_matrix_maximal_dimension_access) {
706 Dim_opt::update_up(dim == Master_matrix::template get_null_value<Dimension>()
707 ? (boundary.size() == 0 ? 0 : boundary.size() - 1)
708 : dim);
709 }
710
711 return _reduce_boundary(cellID, boundary, dim);
712}
713
714template <class Master_matrix>
716{
717 if constexpr (Master_matrix::Option_list::has_map_column_container) {
718 return matrix_.at(columnIndex);
719 } else {
720 return matrix_[columnIndex];
721 }
722}
723
724template <class Master_matrix>
726 Index columnIndex) const
727{
728 if constexpr (Master_matrix::Option_list::has_map_column_container) {
729 return matrix_.at(columnIndex);
730 } else {
731 return matrix_[columnIndex];
732 }
733}
734
735template <class Master_matrix>
737{
738 static_assert(Master_matrix::Option_list::has_removable_columns,
739 "'remove_maximal_cell' is not implemented for the chosen options.");
740 static_assert(Master_matrix::Option_list::has_map_column_container &&
741 Master_matrix::Option_list::has_vine_update &&
742 Master_matrix::Option_list::has_column_pairings,
743 "'remove_maximal_cell' is not implemented for the chosen options.");
744
745 // TODO: find simple test to verify that col at columnIndex is maximal even without row access.
746
747 const auto& pivotToPosition = Swap_opt::CP::pivotToPosition_;
748 auto it = pivotToPosition.find(cellID);
749 if (it == pivotToPosition.end()) return; // cell does not exists. TODO: put an assert instead?
750 Pos_index startPos = it->second;
751 Index startIndex = pivotToColumnIndex_.at(cellID);
752
753 if (startPos != _nextPosition() - 1) {
754 std::vector<Index> colToSwap;
755 colToSwap.reserve(matrix_.size());
756
757 for (auto& p : pivotToPosition) {
758 if (p.second > startPos) colToSwap.push_back(pivotToColumnIndex_.at(p.first));
759 }
760 std::sort(colToSwap.begin(), colToSwap.end(), [&](Index c1, Index c2) {
761 return pivotToPosition.at(get_pivot(c1)) < pivotToPosition.at(get_pivot(c2));
762 });
763
764 for (Index i : colToSwap) {
765 startIndex = Swap_opt::vine_swap(startIndex, i);
766 }
767 }
768
769 _remove_last(startIndex);
770}
771
772template <class Master_matrix>
774 const std::vector<ID_index>& columnsToSwap)
775{
776 static_assert(Master_matrix::Option_list::has_removable_columns,
777 "'remove_maximal_cell' is not implemented for the chosen options.");
778 static_assert(Master_matrix::Option_list::has_map_column_container && Master_matrix::Option_list::has_vine_update,
779 "'remove_maximal_cell' is not implemented for the chosen options.");
780
781 // TODO: find simple test to verify that col at columnIndex is maximal even without row access.
782
783 Index startIndex = pivotToColumnIndex_.at(cellID);
784
785 for (ID_index i : columnsToSwap) {
786 startIndex = Swap_opt::vine_swap(startIndex, pivotToColumnIndex_.at(i));
787 }
788
789 _remove_last(startIndex);
790}
791
792template <class Master_matrix>
794{
795 static_assert(Master_matrix::Option_list::has_removable_columns,
796 "'remove_last' is not implemented for the chosen options.");
797 static_assert(Master_matrix::Option_list::has_map_column_container || !Master_matrix::Option_list::has_vine_update,
798 "'remove_last' is not implemented for the chosen options.");
799
800 if (nextIndex_ == 0 || matrix_.empty()) return; // empty matrix
801
802 if constexpr (Master_matrix::Option_list::has_vine_update) {
803 // careful: linear because of the search of the last index. It is better to keep track of the @ref IDIdx index
804 // of the last column while performing swaps (or the @ref MatIdx with the return values of `vine_swap` + get_pivot)
805 // and then call `remove_maximal_cell` with it and an empty `columnsToSwap`.
806
807 ID_index pivot = 0;
808 Index colIndex = 0;
809 for (auto& p : pivotToColumnIndex_) {
810 if (p.first > pivot) { // pivots have to be strictly increasing in order of filtration
811 pivot = p.first;
812 colIndex = p.second;
813 }
814 }
815 _remove_last(colIndex);
816 } else {
817 _remove_last(nextIndex_ - 1);
818 }
819}
820
821template <class Master_matrix>
823{
824 if constexpr (Master_matrix::Option_list::has_map_column_container) {
825 return matrix_.size();
826 } else {
827 return nextIndex_; // matrix could have been resized much bigger while insert
828 }
829}
830
831template <class Master_matrix>
833 Index columnIndex) const
834{
835 return get_column(columnIndex).get_dimension();
836}
837
838template <class Master_matrix>
839inline void Chain_matrix<Master_matrix>::add_to(Index sourceColumnIndex, Index targetColumnIndex)
840{
841 auto& col = get_column(targetColumnIndex);
842 _add_to(col, [&]() { col += get_column(sourceColumnIndex); });
843}
844
845template <class Master_matrix>
847 const Field_element& coefficient,
848 Index targetColumnIndex)
849{
850 auto& col = get_column(targetColumnIndex);
851 _add_to(col, [&]() { col.multiply_target_and_add(coefficient, get_column(sourceColumnIndex)); });
852}
853
854template <class Master_matrix>
856 Index sourceColumnIndex,
857 Index targetColumnIndex)
858{
859 auto& col = get_column(targetColumnIndex);
860 _add_to(col, [&]() { col.multiply_source_and_add(get_column(sourceColumnIndex), coefficient); });
861}
862
863template <class Master_matrix>
864inline bool Chain_matrix<Master_matrix>::is_zero_entry(Index columnIndex, ID_index rowIndex) const
865{
866 return !get_column(columnIndex).is_non_zero(rowIndex);
867}
868
869template <class Master_matrix>
871{
872 return get_column(columnIndex).is_empty();
873}
874
875template <class Master_matrix>
877 ID_index cellID) const
878{
879 if constexpr (Master_matrix::Option_list::has_map_column_container) {
880 return pivotToColumnIndex_.at(cellID);
881 } else {
882 return pivotToColumnIndex_[cellID];
883 }
884}
885
886template <class Master_matrix>
888{
889 return get_column(columnIndex).get_pivot();
890}
891
892template <class Master_matrix>
894{
895 Dim_opt::operator=(other);
896 Swap_opt::operator=(other);
897 Pair_opt::operator=(other);
898 Rep_opt::operator=(other);
899 matrix_.clear();
900 pivotToColumnIndex_ = other.pivotToColumnIndex_;
901 nextIndex_ = other.nextIndex_;
902 colSettings_ = other.colSettings_;
903
904 matrix_.reserve(other.matrix_.size());
905 for (const auto& cont : other.matrix_){
906 if constexpr (Master_matrix::Option_list::has_map_column_container){
907 _container_insert(cont.second, cont.first);
908 } else {
909 _container_insert(cont);
910 }
911 }
912
913 return *this;
914}
915
916template <class Master_matrix>
917inline void Chain_matrix<Master_matrix>::print() const
918{
919 std::cout << "Column Matrix:\n";
920 if constexpr (!Master_matrix::Option_list::has_map_column_container) {
921 for (ID_index i = 0; i < pivotToColumnIndex_.size(); ++i) {
922 Index pos = pivotToColumnIndex_[i];
923 if (pos != Master_matrix::template get_null_value<Index>()){
924 const Column& col = matrix_[pos];
925 for (const auto& entry : col) {
926 std::cout << entry.get_row_index() << " ";
927 }
928 std::cout << "(" << i << ", " << pos << ")\n";
929 }
930 }
931 if constexpr (Master_matrix::Option_list::has_row_access) {
932 std::cout << "\n";
933 std::cout << "Row Matrix:\n";
934 for (ID_index i = 0; i < pivotToColumnIndex_.size(); ++i) {
935 Index pos = pivotToColumnIndex_[i];
936 if (pos != Master_matrix::template get_null_value<Index>()){
937 const Row& row = RA_opt::get_row(pos);
938 for (const auto& entry : row) {
939 std::cout << entry.get_column_index() << " ";
940 }
941 std::cout << "(" << i << ", " << pos << ")\n";
942 }
943 }
944 }
945 } else {
946 for (const auto& p : pivotToColumnIndex_) {
947 const Column& col = matrix_.at(p.second);
948 for (const auto& entry : col) {
949 std::cout << entry.get_row_index() << " ";
950 }
951 std::cout << "(" << p.first << ", " << p.second << ")\n";
952 }
953 if constexpr (Master_matrix::Option_list::has_row_access) {
954 std::cout << "\n";
955 std::cout << "Row Matrix:\n";
956 for (const auto& p : pivotToColumnIndex_) {
957 const Row& row = RA_opt::get_row(p.first);
958 for (const auto& entry : row) {
959 std::cout << entry.get_column_index() << " ";
960 }
961 std::cout << "(" << p.first << ", " << p.second << ")\n";
962 }
963 }
964 }
965 std::cout << "\n";
966}
967
968template <class Master_matrix>
969template <class Boundary_range>
970inline std::vector<typename Master_matrix::Entry_representative> Chain_matrix<Master_matrix>::_reduce_boundary(
971 ID_index cellID, const Boundary_range& boundary, Dimension dim)
972{
973 Tmp_column column(boundary.begin(), boundary.end());
974 if (dim == Master_matrix::template get_null_value<Dimension>())
975 dim = boundary.begin() == boundary.end() ? 0 : boundary.size() - 1;
976 std::vector<Entry_representative> chainsInH; // for corresponding indices in H (paired columns)
977 std::vector<Entry_representative> chainsInF; // for corresponding indices in F (unpaired, essential columns)
978
979 auto get_last = [&column]() {
980 if constexpr (Master_matrix::Option_list::is_z2)
981 return *(column.rbegin());
982 else
983 return column.rbegin()->first;
984 };
985
986 if (boundary.begin() == boundary.end()) {
987 if constexpr (Master_matrix::Option_list::is_z2)
988 column.insert(cellID);
989 else
990 column.emplace(cellID, 1);
991 _insert_chain(column, dim);
992 return chainsInF;
993 }
994
995 Index currentIndex = get_column_with_pivot(get_last());
996
997 while (get_column(currentIndex).is_paired()) {
998 _reduce_by_G(column, chainsInH, currentIndex);
999
1000 if (column.empty()) {
1001 // produce the sum of all col_h in chains_in_H
1002 _build_from_H(cellID, column, chainsInH);
1003 // create a new cycle (in F) sigma - \sum col_h
1004 _insert_chain(column, dim);
1005 return chainsInF;
1006 }
1007
1008 currentIndex = get_column_with_pivot(get_last());
1009 }
1010
1011 while (!column.empty()) {
1012 currentIndex = get_column_with_pivot(get_last());
1013
1014 if (!get_column(currentIndex).is_paired()) {
1015 // only fills currentEssentialCycleIndices if Z2 coefficients, so chainsInF remains empty
1016 _reduce_by_F(column, chainsInF, currentIndex);
1017 } else {
1018 _reduce_by_G(column, chainsInH, currentIndex);
1019 }
1020 }
1021
1022 _update_largest_death_in_F(chainsInF);
1023
1024 // Compute the new column cellID + \sum col_h, for col_h in chains_in_H
1025 _build_from_H(cellID, column, chainsInH);
1026
1027 // Create and insert (\sum col_h) + sigma (in H, paired with chain_fp) in matrix_
1028 if constexpr (Master_matrix::Option_list::is_z2)
1029 _insert_chain(column, dim, chainsInF[0]);
1030 else
1031 _insert_chain(column, dim, chainsInF[0].first);
1032
1033 return chainsInF;
1034}
1035
1036template <class Master_matrix>
1037inline void Chain_matrix<Master_matrix>::_reduce_by_G(Tmp_column& column,
1038 std::vector<Entry_representative>& chainsInH,
1039 Index currentIndex)
1040{
1041 Column& col = get_column(currentIndex);
1042 if constexpr (Master_matrix::Option_list::is_z2) {
1043 _add_to(col, column, 1u); // Reduce with the column col_g
1044 chainsInH.push_back(col.get_paired_chain_index()); // keep the col_h with which col_g is paired
1045 } else {
1046 Field_element coef = col.get_pivot_value();
1047 auto& operators = colSettings_->operators;
1048 coef = operators.get_inverse(coef);
1049 operators.multiply_inplace(coef, operators.get_characteristic() - column.rbegin()->second);
1050
1051 _add_to(col, column, coef); // Reduce with the column col_g
1052 chainsInH.emplace_back(col.get_paired_chain_index(), coef); // keep the col_h with which col_g is paired
1053 }
1054}
1055
1056template <class Master_matrix>
1057inline void Chain_matrix<Master_matrix>::_reduce_by_F(Tmp_column& column,
1058 std::vector<Entry_representative>& chainsInF,
1059 Index currentIndex)
1060{
1061 Column& col = get_column(currentIndex);
1062 if constexpr (Master_matrix::Option_list::is_z2) {
1063 _add_to(col, column, 1u); // Reduce with the column col_g
1064 chainsInF.push_back(currentIndex);
1065 } else {
1066 Field_element coef = col.get_pivot_value();
1067 auto& operators = colSettings_->operators;
1068 coef = operators.get_inverse(coef);
1069 operators.multiply_inplace(coef, operators.get_characteristic() - column.rbegin()->second);
1070
1071 _add_to(col, column, coef); // Reduce with the column col_g
1072 chainsInF.emplace_back(currentIndex, operators.get_characteristic() - coef);
1073 }
1074}
1075
1076template <class Master_matrix>
1077inline void Chain_matrix<Master_matrix>::_build_from_H(ID_index cellID,
1078 Tmp_column& column,
1079 std::vector<Entry_representative>& chainsInH)
1080{
1081 if constexpr (Master_matrix::Option_list::is_z2) {
1082 column.insert(cellID);
1083 for (Index idx_h : chainsInH) {
1084 _add_to(get_column(idx_h), column, 1u);
1085 }
1086 } else {
1087 column.emplace(cellID, 1);
1088 for (std::pair<Index, Field_element>& idx_h : chainsInH) {
1089 _add_to(get_column(idx_h.first), column, idx_h.second);
1090 }
1091 }
1092}
1093
1094template <class Master_matrix>
1095inline void Chain_matrix<Master_matrix>::_update_largest_death_in_F(const std::vector<Entry_representative>& chainsInF)
1096{
1097 if constexpr (Master_matrix::Option_list::is_z2) {
1098 Index toUpdate = chainsInF[0];
1099 for (auto other_col_it = chainsInF.begin() + 1; other_col_it != chainsInF.end(); ++other_col_it) {
1100 add_to(*other_col_it, toUpdate);
1101 }
1102 } else {
1103 Index toUpdate = chainsInF[0].first;
1104 get_column(toUpdate) *= chainsInF[0].second;
1105 for (auto other_col_it = chainsInF.begin() + 1; other_col_it != chainsInF.end(); ++other_col_it) {
1106 multiply_source_and_add_to(other_col_it->second, other_col_it->first, toUpdate);
1107 }
1108 }
1109}
1110
1111template <class Master_matrix>
1112inline void Chain_matrix<Master_matrix>::_insert_chain(const Tmp_column& column, Dimension dimension)
1113{
1114 _container_insert(column, nextIndex_, dimension);
1115 _add_bar(dimension);
1116
1117 ++nextIndex_;
1118}
1119
1120template <class Master_matrix>
1121inline void Chain_matrix<Master_matrix>::_insert_chain(const Tmp_column& column, Dimension dimension, Index pair)
1122{
1123 // true when no vine updates and if nextIndex_ is updated in remove_last for special case of no vines
1124 // because then @ref PosIdx == @ref MatIdx
1125 Pos_index pairPos = pair;
1126
1127 _container_insert(column, nextIndex_, dimension);
1128
1129 get_column(nextIndex_).assign_paired_chain(pair);
1130 auto& pairCol = get_column(pair);
1131 pairCol.assign_paired_chain(nextIndex_);
1132
1133 if constexpr (Master_matrix::Option_list::has_column_pairings && Master_matrix::Option_list::has_vine_update) {
1134 pairPos = Swap_opt::CP::pivotToPosition_[pairCol.get_pivot()];
1135 }
1136
1137 _update_barcode(pairPos);
1138
1139 ++nextIndex_;
1140}
1141
1142template <class Master_matrix>
1143inline void Chain_matrix<Master_matrix>::_add_to(const Column& column,
1144 Tmp_column& set,
1145 [[maybe_unused]] unsigned int coef)
1146{
1147 if constexpr (Master_matrix::Option_list::is_z2) {
1148 std::pair<typename std::set<Index>::iterator, bool> res_insert;
1149 for (const Entry& entry : column) {
1150 res_insert = set.insert(entry.get_row_index());
1151 if (!res_insert.second) {
1152 set.erase(res_insert.first);
1153 }
1154 }
1155 } else {
1156 auto& operators = colSettings_->operators;
1157 for (const Entry& entry : column) {
1158 auto res = set.emplace(entry.get_row_index(), entry.get_element());
1159 if (res.second){
1160 operators.multiply_inplace(res.first->second, coef);
1161 } else {
1162 operators.multiply_and_add_inplace_back(entry.get_element(), coef, res.first->second);
1163 if (res.first->second == Field_operators::get_additive_identity()) {
1164 set.erase(res.first);
1165 }
1166 }
1167 }
1168 }
1169}
1170
1171template <class Master_matrix>
1172template <typename F>
1173inline void Chain_matrix<Master_matrix>::_add_to(Column& target, F&& addition)
1174{
1175 auto pivot = target.get_pivot();
1176 addition();
1177
1178 if (pivot != target.get_pivot()) {
1179 if constexpr (Master_matrix::Option_list::has_map_column_container) {
1180 std::swap(pivotToColumnIndex_.at(pivot), pivotToColumnIndex_.at(target.get_pivot()));
1181 } else {
1182 std::swap(pivotToColumnIndex_[pivot], pivotToColumnIndex_[target.get_pivot()]);
1183 }
1184 }
1185}
1186
1187template <class Master_matrix>
1188inline void Chain_matrix<Master_matrix>::_remove_last(Index lastIndex)
1189{
1190 static_assert(Master_matrix::Option_list::has_removable_columns,
1191 "'_remove_last' is not implemented for the chosen options.");
1192 static_assert(Master_matrix::Option_list::has_map_column_container || !Master_matrix::Option_list::has_vine_update,
1193 "'_remove_last' is not implemented for the chosen options.");
1194
1195 ID_index pivot;
1196
1197 if constexpr (Master_matrix::Option_list::has_map_column_container) {
1198 auto itToErase = matrix_.find(lastIndex);
1199 Column& colToErase = itToErase->second;
1200 pivot = colToErase.get_pivot();
1201
1202 if constexpr (Master_matrix::Option_list::has_matrix_maximal_dimension_access) {
1203 Dim_opt::update_down(colToErase.get_dimension());
1204 }
1205
1206 if (colToErase.is_paired()) matrix_.at(colToErase.get_paired_chain_index()).unassign_paired_chain();
1207 pivotToColumnIndex_.erase(pivot);
1208 matrix_.erase(itToErase);
1209 } else {
1210 GUDHI_CHECK(lastIndex == nextIndex_ - 1 && nextIndex_ == matrix_.size(),
1211 std::logic_error("Chain_matrix::_remove_last - Indexation problem."));
1212
1213 Column& colToErase = matrix_[lastIndex];
1214 pivot = colToErase.get_pivot();
1215
1216 if constexpr (Master_matrix::Option_list::has_matrix_maximal_dimension_access) {
1217 Dim_opt::update_down(colToErase.get_dimension());
1218 }
1219
1220 if (colToErase.is_paired()) matrix_.at(colToErase.get_paired_chain_index()).unassign_paired_chain();
1221 pivotToColumnIndex_[pivot] = Master_matrix::template get_null_value<Index>();
1222 matrix_.pop_back();
1223 // TODO: resize matrix_ when a lot is removed? Could be not the best strategy if user inserts a lot back afterwards.
1224 }
1225
1226 if constexpr (!Master_matrix::Option_list::has_vine_update) {
1227 --nextIndex_; // should not be updated when there are vine updates, as possibly lastIndex != nextIndex - 1
1228 }
1229
1230 if constexpr (Master_matrix::Option_list::has_column_pairings) {
1231 auto it = _indexToBar().find(--_nextPosition());
1232 typename Barcode::iterator bar = it->second;
1233
1234 if (bar->death == Master_matrix::template get_null_value<Pos_index>())
1235 _barcode().erase(bar);
1236 else
1237 bar->death = Master_matrix::template get_null_value<Pos_index>();
1238
1239 _indexToBar().erase(it);
1240 if constexpr (Master_matrix::Option_list::has_vine_update) Swap_opt::CP::pivotToPosition_.erase(pivot);
1241 }
1242
1243 if constexpr (Master_matrix::Option_list::has_row_access) {
1244 GUDHI_CHECK(
1245 RA_opt::get_row(pivot).size() == 0,
1246 std::invalid_argument(
1247 "Chain_matrix::_remove_last - Column asked to be removed does not corresponds to a maximal simplex."));
1248 if constexpr (Master_matrix::Option_list::has_removable_rows) {
1249 RA_opt::erase_empty_row(pivot);
1250 }
1251 }
1252}
1253
1254template <class Master_matrix>
1255inline void Chain_matrix<Master_matrix>::_update_barcode(Pos_index birth)
1256{
1257 if constexpr (Master_matrix::Option_list::has_column_pairings) {
1258 if constexpr (Master_matrix::Option_list::has_removable_columns) {
1259 auto& barIt = _indexToBar().at(birth);
1260 barIt->death = _nextPosition();
1261 _indexToBar().try_emplace(_nextPosition(), barIt); // list so iterators are stable
1262 } else {
1263 _barcode()[_indexToBar()[birth]].death = _nextPosition();
1264 _indexToBar().push_back(_indexToBar()[birth]);
1265 }
1266 ++_nextPosition();
1267 }
1268}
1269
1270template <class Master_matrix>
1271inline void Chain_matrix<Master_matrix>::_add_bar(Dimension dim)
1272{
1273 if constexpr (Master_matrix::Option_list::has_column_pairings) {
1274 _barcode().emplace_back(_nextPosition(), Master_matrix::template get_null_value<Pos_index>(), dim);
1275 if constexpr (Master_matrix::Option_list::has_removable_columns) {
1276 _indexToBar().try_emplace(_nextPosition(), --_barcode().end());
1277 } else {
1278 _indexToBar().push_back(_barcode().size() - 1);
1279 }
1280 ++_nextPosition();
1281 }
1282}
1283
1284template <class Master_matrix>
1285template <class Container>
1286inline void Chain_matrix<Master_matrix>::_container_insert(const Container& column, Index pos, Dimension dim)
1287{
1288 ID_index pivot;
1289 if constexpr (Master_matrix::Option_list::is_z2) {
1290 pivot = *(column.rbegin());
1291 } else {
1292 pivot = column.rbegin()->first;
1293 }
1294 if constexpr (Master_matrix::Option_list::has_map_column_container) {
1295 pivotToColumnIndex_.try_emplace(pivot, pos);
1296 if constexpr (Master_matrix::Option_list::has_row_access) {
1297 matrix_.try_emplace(pos, Column(pos, column, dim, RA_opt::rows_, colSettings_));
1298 } else {
1299 matrix_.try_emplace(pos, Column(column, dim, colSettings_));
1300 }
1301 } else {
1302 if constexpr (Master_matrix::Option_list::has_row_access) {
1303 matrix_.emplace_back(pos, column, dim, RA_opt::rows_, colSettings_);
1304 } else {
1305 matrix_.emplace_back(column, dim, colSettings_);
1306 }
1307 pivotToColumnIndex_[pivot] = pos;
1308 }
1309}
1310
1311template <class Master_matrix>
1312inline void Chain_matrix<Master_matrix>::_container_insert(const Column& column, [[maybe_unused]] Index pos)
1313{
1314 if constexpr (Master_matrix::Option_list::has_map_column_container) {
1315 if constexpr (Master_matrix::Option_list::has_row_access) {
1316 matrix_.try_emplace(pos, Column(column, column.get_column_index(), RA_opt::rows_, colSettings_));
1317 } else {
1318 matrix_.try_emplace(pos, Column(column, colSettings_));
1319 }
1320 } else {
1321 if constexpr (Master_matrix::Option_list::has_row_access) {
1322 matrix_.emplace_back(column, column.get_column_index(), RA_opt::rows_, colSettings_);
1323 } else {
1324 matrix_.emplace_back(column, colSettings_);
1325 }
1326 }
1327}
1328
1329template <class Master_matrix>
1331{
1332 if constexpr (Master_matrix::Option_list::has_vine_update)
1333 return Swap_opt::template Chain_pairing<Master_matrix>::barcode_;
1334 else
1335 return Pair_opt::barcode_;
1336}
1337
1338template <class Master_matrix>
1339inline constexpr typename Chain_matrix<Master_matrix>::Bar_dictionary&
1341{
1342 if constexpr (Master_matrix::Option_list::has_vine_update)
1343 return Swap_opt::template Chain_pairing<Master_matrix>::indexToBar_;
1344 else
1345 return Pair_opt::indexToBar_;
1346}
1347
1348template <class Master_matrix>
1350{
1351 if constexpr (Master_matrix::Option_list::has_vine_update)
1352 return Swap_opt::template Chain_pairing<Master_matrix>::nextPosition_;
1353 else
1354 return Pair_opt::nextPosition_;
1355}
1356
1357} // namespace persistence_matrix
1358} // namespace Gudhi
1359
1360#endif // PM_CHAIN_MATRIX_H
Contains the Gudhi::persistence_matrix::Id_to_index_overlay class.
Matrix structure storing a compatible base of a filtered chain complex. See . The base is constructed...
Definition: Chain_matrix.h:50
void remove_maximal_cell(ID_index cellID)
Only available if PersistenceMatrixOptions::has_removable_columns and PersistenceMatrixOptions::has_v...
Definition: Chain_matrix.h:736
friend void swap(Chain_matrix &matrix1, Chain_matrix &matrix2)
Swap operator.
Definition: Chain_matrix.h:459
std::vector< Entry_representative > insert_boundary(const Boundary_range &boundary, Dimension dim=Master_matrix::template get_null_value< Dimension >())
Inserts at the end of the matrix a new ordered column corresponding to the given boundary....
Column & get_column(Index columnIndex)
Returns the column at the given MatIdx index. The type of the column depends on the chosen options,...
Definition: Chain_matrix.h:715
void reset(Column_settings *colSettings)
Resets the matrix to an empty matrix.
Definition: Chain_matrix.h:445
typename Master_matrix::Index Index
Definition: Chain_matrix.h:66
void multiply_source_and_add_to(const Field_element &coefficient, Index sourceColumnIndex, Index targetColumnIndex)
Multiplies the source column with the coefficient before adding it to the target column....
Definition: Chain_matrix.h:855
typename Master_matrix::Pos_index Pos_index
Definition: Chain_matrix.h:68
bool is_zero_entry(Index columnIndex, ID_index rowIndex) const
Indicates if the entry at given coordinates has value zero.
Definition: Chain_matrix.h:864
typename Master_matrix::Entry_constructor Entry_constructor
Definition: Chain_matrix.h:61
void multiply_target_and_add_to(Index sourceColumnIndex, const Field_element &coefficient, Index targetColumnIndex)
Multiplies the target column with the coefficient and then adds the source column to it....
Definition: Chain_matrix.h:846
Index get_column_with_pivot(ID_index cellID) const
Returns the column with given row index as pivot. Assumes that the pivot exists.
Definition: Chain_matrix.h:876
Chain_matrix(Column_settings *colSettings)
Constructs an empty matrix. Only available if PersistenceMatrixOptions::has_column_pairings is true o...
Definition: Chain_matrix.h:527
typename Master_matrix::Row Row
Definition: Chain_matrix.h:59
typename Master_matrix::Field_operators Field_operators
Field operators class. Necessary only if PersistenceMatrixOptions::is_z2 is false.
Definition: Chain_matrix.h:55
typename Master_matrix::Element Field_element
Definition: Chain_matrix.h:56
void add_to(Index sourceColumnIndex, Index targetColumnIndex)
Adds column at sourceColumnIndex onto the column at targetColumnIndex in the matrix.
Definition: Chain_matrix.h:839
ID_index get_pivot(Index columnIndex)
Returns the row index of the pivot of the given column.
Definition: Chain_matrix.h:887
typename Master_matrix::Entry_representative Entry_representative
Definition: Chain_matrix.h:65
bool is_zero_column(Index columnIndex)
Indicates if the column at given index has value zero. Note that if the matrix is valid,...
Definition: Chain_matrix.h:870
Index get_number_of_columns() const
Returns the current number of columns in the matrix.
Definition: Chain_matrix.h:822
std::vector< Entry_representative > insert_boundary(ID_index cellID, const Boundary_range &boundary, Dimension dim=Master_matrix::template get_null_value< Dimension >())
It does the same as the other version, but allows the boundary cells to be identified without restric...
Dimension get_column_dimension(Index columnIndex) const
Returns the dimension of the given column.
Definition: Chain_matrix.h:832
void remove_last()
Only available if PersistenceMatrixOptions::has_removable_columns is true and, if PersistenceMatrixOp...
Definition: Chain_matrix.h:793
typename Master_matrix::ID_index ID_index
Definition: Chain_matrix.h:67
typename Master_matrix::Column Column
Definition: Chain_matrix.h:57
typename Master_matrix::Matrix_entry Entry
Definition: Chain_matrix.h:60
typename Master_matrix::Dimension Dimension
Definition: Chain_matrix.h:69
Chain_matrix & operator=(const Chain_matrix &other)
Assign operator.
Definition: Chain_matrix.h:893
typename Master_matrix::Boundary Boundary
Definition: Chain_matrix.h:64
typename Master_matrix::Column_settings Column_settings
Definition: Chain_matrix.h:63
Overlay for non-basic matrices replacing all input and output MatIdx indices of the original methods ...
Definition: Id_to_index_overlay.h:42
Data structure for matrices, and in particular thought for matrices representing filtered complexes i...
Definition: Matrix.h:144
typename std::conditional< hasFixedBarcode, std::vector< Bar >, typename std::conditional< PersistenceMatrixOptions::has_removable_columns, std::list< Bar >, std::vector< Bar > >::type >::type Barcode
Type of the computed barcode. It is either a list of Matrix::Bar or a vector of Matrix::Bar,...
Definition: Matrix.h:447
ID_index get_pivot(Index columnIndex)
Returns the row index of the pivot of the given column. Only available for non-basic matrices.
Definition: Matrix.h:1930
typename PersistenceMatrixOptions::Dimension Dimension
Definition: Matrix.h:150
typename PersistenceMatrixOptions::Index Index
Definition: Matrix.h:147
typename PersistenceMatrixOptions::Index ID_index
Definition: Matrix.h:148
typename PersistenceMatrixOptions::Index Pos_index
Definition: Matrix.h:149
typename std::conditional< PersistenceMatrixOptions::column_type==Column_types::HEAP, Matrix_heap_column, typename std::conditional< PersistenceMatrixOptions::column_type==Column_types::LIST, Matrix_list_column, typename std::conditional< PersistenceMatrixOptions::column_type==Column_types::SET, Matrix_set_column, typename std::conditional< PersistenceMatrixOptions::column_type==Column_types::UNORDERED_SET, Matrix_unordered_set_column, typename std::conditional< PersistenceMatrixOptions::column_type==Column_types::VECTOR, Matrix_vector_column, typename std::conditional< PersistenceMatrixOptions::column_type==Column_types::INTRUSIVE_LIST, Matrix_intrusive_list_column, typename std::conditional< PersistenceMatrixOptions::column_type==Column_types::NAIVE_VECTOR, Matrix_naive_vector_column, typename std::conditional< PersistenceMatrixOptions::column_type==Column_types::SMALL_VECTOR, Matrix_small_vector_column, Matrix_intrusive_set_column >::type >::type >::type >::type >::type >::type >::type >::type Column
Type of the columns stored in the matrix. The type depends on the value of PersistenceMatrixOptions::...
Definition: Matrix.h:368
Gudhi namespace.
Definition: SimplicialComplexForAlpha.h:14
static const bool has_map_column_container
If set to true, the underlying container containing the matrix columns is an std::unordered_map....
Definition: PersistenceMatrixOptions.h:101
static const bool has_row_access
If set to true, enables the method Matrix::get_row.
Definition: PersistenceMatrixOptions.h:111
static const bool is_z2
If true, indicates that the values contained in the matrix are in and can therefore be treated like ...
Definition: PersistenceMatrixOptions.h:56