Loading...
Searching...
No Matches
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_type = typename Master_matrix::element_type;
57 using Column_type = typename Master_matrix::Column_type;
58 using Row_type = typename Master_matrix::Row_type;
60 using Cell = typename Master_matrix::Cell_type;
61 using Cell_constructor = typename Master_matrix::Cell_constructor;
62 using Column_settings = typename Master_matrix::Column_settings;
64 using boundary_type = typename Master_matrix::boundary_type;
65 using cell_rep_type = typename Master_matrix::cell_rep_type;
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_type = typename Master_matrix::dimension_type;
79 Chain_matrix(Column_settings* colSettings);
104 template <class Boundary_type = boundary_type>
105 Chain_matrix(const std::vector<Boundary_type>& 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);
177 template <typename BirthComparatorFunction, typename DeathComparatorFunction, class Boundary_type = boundary_type>
178 Chain_matrix(const std::vector<Boundary_type>& orderedBoundaries,
179 Column_settings* colSettings,
180 const BirthComparatorFunction& birthComparator,
181 const DeathComparatorFunction& deathComparator);
204 template <typename BirthComparatorFunction, typename DeathComparatorFunction>
205 Chain_matrix(unsigned int numberOfColumns,
206 Column_settings* colSettings,
207 const BirthComparatorFunction& birthComparator,
208 const DeathComparatorFunction& deathComparator);
218 Chain_matrix(const Chain_matrix& matrixToCopy,
219 Column_settings* colSettings = nullptr);
225 Chain_matrix(Chain_matrix&& other) noexcept;
226
247 template <class Boundary_type = boundary_type>
248 std::vector<cell_rep_type> insert_boundary(const Boundary_type& boundary, dimension_type dim = -1);
266 template <class Boundary_type = boundary_type>
267 std::vector<cell_rep_type> insert_boundary(id_index faceID, const Boundary_type& boundary, dimension_type dim = -1);
275 Column_type& get_column(index columnIndex);
283 const Column_type& get_column(index columnIndex) const;
300 void remove_maximal_face(id_index faceID);
324 void remove_maximal_face(id_index faceID, const std::vector<id_index>& columnsToSwap);
339 void remove_last();
340
347
354 dimension_type get_column_dimension(index columnIndex) const;
355
366 void add_to(index sourceColumnIndex, index targetColumnIndex);
379 void multiply_target_and_add_to(index sourceColumnIndex,
380 const Field_element_type& coefficient,
381 index targetColumnIndex);
394 void multiply_source_and_add_to(const Field_element_type& coefficient,
395 index sourceColumnIndex,
396 index targetColumnIndex);
397
406 bool is_zero_cell(index columnIndex, id_index rowIndex) const;
415 bool is_zero_column(index columnIndex);
416
430 id_index get_pivot(index columnIndex);
431
438 void reset(Column_settings* colSettings) {
439 matrix_.clear();
440 pivotToColumnIndex_.clear();
441 nextIndex_ = 0;
442 colSettings_ = colSettings;
443 }
444
448 Chain_matrix& operator=(const Chain_matrix& other);
452 friend void swap(Chain_matrix& matrix1, Chain_matrix& matrix2) {
453 swap(static_cast<typename Master_matrix::Matrix_dimension_option&>(matrix1),
454 static_cast<typename Master_matrix::Matrix_dimension_option&>(matrix2));
455 swap(static_cast<typename Master_matrix::Chain_pairing_option&>(matrix1),
456 static_cast<typename Master_matrix::Chain_pairing_option&>(matrix2));
457 swap(static_cast<typename Master_matrix::Chain_vine_swap_option&>(matrix1),
458 static_cast<typename Master_matrix::Chain_vine_swap_option&>(matrix2));
459 swap(static_cast<typename Master_matrix::Chain_representative_cycles_option&>(matrix1),
460 static_cast<typename Master_matrix::Chain_representative_cycles_option&>(matrix2));
461 matrix1.matrix_.swap(matrix2.matrix_);
462 matrix1.pivotToColumnIndex_.swap(matrix2.pivotToColumnIndex_);
463 std::swap(matrix1.nextIndex_, matrix2.nextIndex_);
464 std::swap(matrix1.colSettings_, matrix2.colSettings_);
465
466 if constexpr (Master_matrix::Option_list::has_row_access) {
467 swap(static_cast<typename Master_matrix::Matrix_row_access_option&>(matrix1),
468 static_cast<typename Master_matrix::Matrix_row_access_option&>(matrix2));
469 }
470 }
471
472 void print() const; // for debug
473
474 friend class Id_to_index_overlay<Chain_matrix<Master_matrix>, Master_matrix>;
475
476 private:
477 using dim_opt = typename Master_matrix::Matrix_dimension_option;
478 using swap_opt = typename Master_matrix::Chain_vine_swap_option;
479 using pair_opt = typename Master_matrix::Chain_pairing_option;
480 using rep_opt = typename Master_matrix::Chain_representative_cycles_option;
481 using ra_opt = typename Master_matrix::Matrix_row_access_option;
482 using matrix_type = typename Master_matrix::column_container_type;
483 using dictionnary_type = typename Master_matrix::template dictionnary_type<index>;
484 using barcode_type = typename Master_matrix::barcode_type;
485 using bar_dictionnary_type = typename Master_matrix::bar_dictionnary_type;
486 using tmp_column_type = typename std::conditional<
487 Master_matrix::Option_list::is_z2,
488 std::set<id_index>,
489 std::map<id_index, Field_element_type>
490 >::type;
491
492 matrix_type matrix_;
493 dictionnary_type pivotToColumnIndex_;
494 index nextIndex_;
495 Column_settings* colSettings_;
497 template <class Boundary_type>
498 std::vector<cell_rep_type> _reduce_boundary(id_index faceID, const Boundary_type& boundary, dimension_type dim);
499 void _reduce_by_G(tmp_column_type& column, std::vector<cell_rep_type>& chainsInH, index currentPivot);
500 void _reduce_by_F(tmp_column_type& column, std::vector<cell_rep_type>& chainsInF, index currentPivot);
501 void _build_from_H(id_index faceID, tmp_column_type& column, std::vector<cell_rep_type>& chainsInH);
502 void _update_largest_death_in_F(const std::vector<cell_rep_type>& chainsInF);
503 void _insert_chain(const tmp_column_type& column, dimension_type dimension);
504 void _insert_chain(const tmp_column_type& column, dimension_type dimension, index pair);
505 void _add_to(const Column_type& column, tmp_column_type& set, unsigned int coef);
506 template <typename F>
507 void _add_to(Column_type& target, F&& addition);
508 void _remove_last(index lastIndex);
509 void _update_barcode(pos_index birth);
510 void _add_bar(dimension_type dim);
511 template <class Container_type>
512 void _container_insert(const Container_type& column, index pos, dimension_type dim);
513 void _container_insert(const Column_type& column, [[maybe_unused]] index pos = 0);
514
515 constexpr barcode_type& _barcode();
516 constexpr bar_dictionnary_type& _indexToBar();
517 constexpr pos_index& _nextPosition();
518};
519
520template <class Master_matrix>
522 : dim_opt(-1),
523 pair_opt(),
524 swap_opt(),
525 rep_opt(),
526 ra_opt(),
527 nextIndex_(0),
528 colSettings_(colSettings)
529{}
530
531template <class Master_matrix>
532template <class Boundary_type>
533inline Chain_matrix<Master_matrix>::Chain_matrix(const std::vector<Boundary_type>& orderedBoundaries,
534 Column_settings* colSettings)
535 : dim_opt(-1),
536 pair_opt(),
537 swap_opt(),
538 rep_opt(),
539 ra_opt(orderedBoundaries.size()),
540 nextIndex_(0),
541 colSettings_(colSettings)
542{
543 matrix_.reserve(orderedBoundaries.size());
544 if constexpr (Master_matrix::Option_list::has_map_column_container) {
545 pivotToColumnIndex_.reserve(orderedBoundaries.size());
546 } else {
547 pivotToColumnIndex_.resize(orderedBoundaries.size(), -1);
548 }
549
550 for (const Boundary_type& b : orderedBoundaries) {
552 }
553}
554
555template <class Master_matrix>
556inline Chain_matrix<Master_matrix>::Chain_matrix(unsigned int numberOfColumns,
557 Column_settings* colSettings)
558 : dim_opt(-1),
559 pair_opt(),
560 swap_opt(),
561 rep_opt(),
562 ra_opt(numberOfColumns),
563 nextIndex_(0),
564 colSettings_(colSettings)
565{
566 matrix_.reserve(numberOfColumns);
567 if constexpr (Master_matrix::Option_list::has_map_column_container) {
568 pivotToColumnIndex_.reserve(numberOfColumns);
569 } else {
570 pivotToColumnIndex_.resize(numberOfColumns, -1);
571 }
572}
573
574template <class Master_matrix>
575template <typename BirthComparatorFunction, typename DeathComparatorFunction>
577 const BirthComparatorFunction& birthComparator,
578 const DeathComparatorFunction& deathComparator)
579 : dim_opt(-1),
580 pair_opt(),
581 swap_opt(birthComparator, deathComparator),
582 rep_opt(),
583 ra_opt(),
584 nextIndex_(0),
585 colSettings_(colSettings)
586{}
587
588template <class Master_matrix>
589template <typename BirthComparatorFunction, typename DeathComparatorFunction, class Boundary_type>
590inline Chain_matrix<Master_matrix>::Chain_matrix(const std::vector<Boundary_type>& orderedBoundaries,
591 Column_settings* colSettings,
592 const BirthComparatorFunction& birthComparator,
593 const DeathComparatorFunction& deathComparator)
594 : dim_opt(-1),
595 pair_opt(),
596 swap_opt(birthComparator, deathComparator),
597 rep_opt(),
598 ra_opt(orderedBoundaries.size()),
599 nextIndex_(0),
600 colSettings_(colSettings)
601{
602 matrix_.reserve(orderedBoundaries.size());
603 if constexpr (Master_matrix::Option_list::has_map_column_container) {
604 pivotToColumnIndex_.reserve(orderedBoundaries.size());
605 } else {
606 pivotToColumnIndex_.resize(orderedBoundaries.size(), -1);
607 }
608 for (const Boundary_type& b : orderedBoundaries) {
610 }
611}
612
613template <class Master_matrix>
614template <typename BirthComparatorFunction, typename DeathComparatorFunction>
615inline Chain_matrix<Master_matrix>::Chain_matrix(unsigned int numberOfColumns,
616 Column_settings* colSettings,
617 const BirthComparatorFunction& birthComparator,
618 const DeathComparatorFunction& deathComparator)
619 : dim_opt(-1),
620 pair_opt(),
621 swap_opt(birthComparator, deathComparator),
622 rep_opt(),
623 ra_opt(numberOfColumns),
624 nextIndex_(0),
625 colSettings_(colSettings)
626{
627 matrix_.reserve(numberOfColumns);
628 if constexpr (Master_matrix::Option_list::has_map_column_container) {
629 pivotToColumnIndex_.reserve(numberOfColumns);
630 } else {
631 pivotToColumnIndex_.resize(numberOfColumns, -1);
632 }
633}
634
635template <class Master_matrix>
637 Column_settings* colSettings)
638 : dim_opt(static_cast<const dim_opt&>(matrixToCopy)),
639 pair_opt(static_cast<const pair_opt&>(matrixToCopy)),
640 swap_opt(static_cast<const swap_opt&>(matrixToCopy)),
641 rep_opt(static_cast<const rep_opt&>(matrixToCopy)),
642 ra_opt(static_cast<const ra_opt&>(matrixToCopy)),
643 pivotToColumnIndex_(matrixToCopy.pivotToColumnIndex_),
644 nextIndex_(matrixToCopy.nextIndex_),
645 colSettings_(colSettings == nullptr ? matrixToCopy.colSettings_ : colSettings)
646{
647 matrix_.reserve(matrixToCopy.matrix_.size());
648 for (const auto& cont : matrixToCopy.matrix_){
649 if constexpr (Master_matrix::Option_list::has_map_column_container){
650 _container_insert(cont.second, cont.first);
651 } else {
652 _container_insert(cont);
653 }
654 }
655}
656
657template <class Master_matrix>
659 : dim_opt(std::move(static_cast<dim_opt&>(other))),
660 pair_opt(std::move(static_cast<pair_opt&>(other))),
661 swap_opt(std::move(static_cast<swap_opt&>(other))),
662 rep_opt(std::move(static_cast<rep_opt&>(other))),
663 ra_opt(std::move(static_cast<ra_opt&>(other))),
664 matrix_(std::move(other.matrix_)),
665 pivotToColumnIndex_(std::move(other.pivotToColumnIndex_)),
666 nextIndex_(std::exchange(other.nextIndex_, 0)),
667 colSettings_(std::exchange(other.colSettings_, nullptr))
668{}
669
670template <class Master_matrix>
671template <class Boundary_type>
672inline std::vector<typename Master_matrix::cell_rep_type> Chain_matrix<Master_matrix>::insert_boundary(
673 const Boundary_type& boundary, dimension_type dim)
674{
675 return insert_boundary(nextIndex_, boundary, dim);
676}
677
678template <class Master_matrix>
679template <class Boundary_type>
680inline std::vector<typename Master_matrix::cell_rep_type> Chain_matrix<Master_matrix>::insert_boundary(
681 id_index faceID, const Boundary_type& boundary, dimension_type dim)
682{
683 if constexpr (!Master_matrix::Option_list::has_map_column_container) {
684 if (pivotToColumnIndex_.size() <= faceID) {
685 pivotToColumnIndex_.resize(faceID * 2 + 1, -1);
686 }
687 }
688
689 if constexpr (Master_matrix::Option_list::has_vine_update && Master_matrix::Option_list::has_column_pairings) {
690 if constexpr (Master_matrix::Option_list::has_map_column_container) {
691 swap_opt::CP::pivotToPosition_.try_emplace(faceID, _nextPosition());
692 } else {
693 if (swap_opt::CP::pivotToPosition_.size() <= faceID)
694 swap_opt::CP::pivotToPosition_.resize(pivotToColumnIndex_.size(), -1);
695 swap_opt::CP::pivotToPosition_[faceID] = _nextPosition();
696 }
697 }
698
699 if constexpr (Master_matrix::Option_list::has_matrix_maximal_dimension_access) {
700 dim_opt::update_up(dim == static_cast<dimension_type>(-1) ? (boundary.size() == 0 ? 0 : boundary.size() - 1) : dim);
701 }
702
703 return _reduce_boundary(faceID, boundary, dim);
704}
705
706template <class Master_matrix>
708{
709 if constexpr (Master_matrix::Option_list::has_map_column_container) {
710 return matrix_.at(columnIndex);
711 } else {
712 return matrix_[columnIndex];
713 }
714}
715
716template <class Master_matrix>
718 index columnIndex) const
719{
720 if constexpr (Master_matrix::Option_list::has_map_column_container) {
721 return matrix_.at(columnIndex);
722 } else {
723 return matrix_[columnIndex];
724 }
725}
726
727template <class Master_matrix>
729{
730 static_assert(Master_matrix::Option_list::has_removable_columns,
731 "'remove_maximal_face' is not implemented for the chosen options.");
732 static_assert(Master_matrix::Option_list::has_map_column_container &&
733 Master_matrix::Option_list::has_vine_update &&
734 Master_matrix::Option_list::has_column_pairings,
735 "'remove_maximal_face' is not implemented for the chosen options.");
736
737 // TODO: find simple test to verify that col at columnIndex is maximal even without row access.
738
739 const auto& pivotToPosition = swap_opt::CP::pivotToPosition_;
740 auto it = pivotToPosition.find(faceID);
741 if (it == pivotToPosition.end()) return; // face does not exists. TODO: put an assert instead?
742 pos_index startPos = it->second;
743 index startIndex = pivotToColumnIndex_.at(faceID);
744
745 if (startPos != _nextPosition() - 1) {
746 std::vector<index> colToSwap;
747 colToSwap.reserve(matrix_.size());
748
749 for (auto& p : pivotToPosition) {
750 if (p.second > startPos) colToSwap.push_back(pivotToColumnIndex_.at(p.first));
751 }
752 std::sort(colToSwap.begin(), colToSwap.end(), [&](index c1, index c2) {
753 return pivotToPosition.at(get_pivot(c1)) < pivotToPosition.at(get_pivot(c2));
754 });
755
756 for (index i : colToSwap) {
757 startIndex = swap_opt::vine_swap(startIndex, i);
758 }
759 }
760
761 _remove_last(startIndex);
762}
763
764template <class Master_matrix>
766 const std::vector<id_index>& columnsToSwap)
767{
768 static_assert(Master_matrix::Option_list::has_removable_columns,
769 "'remove_maximal_face' is not implemented for the chosen options.");
770 static_assert(Master_matrix::Option_list::has_map_column_container && Master_matrix::Option_list::has_vine_update,
771 "'remove_maximal_face' is not implemented for the chosen options.");
772
773 // TODO: find simple test to verify that col at columnIndex is maximal even without row access.
774
775 index startIndex = pivotToColumnIndex_.at(faceID);
776
777 for (id_index i : columnsToSwap) {
778 startIndex = swap_opt::vine_swap(startIndex, pivotToColumnIndex_.at(i));
779 }
780
781 _remove_last(startIndex);
782}
783
784template <class Master_matrix>
786{
787 static_assert(Master_matrix::Option_list::has_removable_columns,
788 "'remove_last' is not implemented for the chosen options.");
789 static_assert(Master_matrix::Option_list::has_map_column_container || !Master_matrix::Option_list::has_vine_update,
790 "'remove_last' is not implemented for the chosen options.");
791
792 if (nextIndex_ == 0 || matrix_.empty()) return; // empty matrix
793
794 if constexpr (Master_matrix::Option_list::has_vine_update) {
795 // carefull: linear because of the search of the last index. It is better to keep track of the @ref IDIdx index
796 // of the last column while performing swaps (or the @ref MatIdx with the return values of `vine_swap` + get_pivot)
797 // and then call `remove_maximal_face` with it and an empty `columnsToSwap`.
798
799 id_index pivot = 0;
800 index colIndex = 0;
801 for (auto& p : pivotToColumnIndex_) {
802 if (p.first > pivot) { // pivots have to be strictly increasing in order of filtration
803 pivot = p.first;
804 colIndex = p.second;
805 }
806 }
807 _remove_last(colIndex);
808 } else {
809 _remove_last(nextIndex_ - 1);
810 }
811}
812
813template <class Master_matrix>
815{
816 if constexpr (Master_matrix::Option_list::has_map_column_container) {
817 return matrix_.size();
818 } else {
819 return nextIndex_; // matrix could have been resized much bigger while insert
820 }
821}
822
823template <class Master_matrix>
825 index columnIndex) const
826{
827 return get_column(columnIndex).get_dimension();
828}
829
830template <class Master_matrix>
831inline void Chain_matrix<Master_matrix>::add_to(index sourceColumnIndex, index targetColumnIndex)
832{
833 auto& col = get_column(targetColumnIndex);
834 _add_to(col, [&]() { col += get_column(sourceColumnIndex); });
835}
836
837template <class Master_matrix>
839 const Field_element_type& coefficient,
840 index targetColumnIndex)
841{
842 auto& col = get_column(targetColumnIndex);
843 _add_to(col, [&]() { col.multiply_target_and_add(coefficient, get_column(sourceColumnIndex)); });
844}
845
846template <class Master_matrix>
848 index sourceColumnIndex,
849 index targetColumnIndex)
850{
851 auto& col = get_column(targetColumnIndex);
852 _add_to(col, [&]() { col.multiply_source_and_add(get_column(sourceColumnIndex), coefficient); });
853}
854
855template <class Master_matrix>
856inline bool Chain_matrix<Master_matrix>::is_zero_cell(index columnIndex, id_index rowIndex) const
857{
858 return !get_column(columnIndex).is_non_zero(rowIndex);
859}
860
861template <class Master_matrix>
863{
864 return get_column(columnIndex).is_empty();
865}
866
867template <class Master_matrix>
869 id_index faceID) const
870{
871 if constexpr (Master_matrix::Option_list::has_map_column_container) {
872 return pivotToColumnIndex_.at(faceID);
873 } else {
874 return pivotToColumnIndex_[faceID];
875 }
876}
877
878template <class Master_matrix>
880{
881 return get_column(columnIndex).get_pivot();
882}
883
884template <class Master_matrix>
886{
887 dim_opt::operator=(other);
888 swap_opt::operator=(other);
889 pair_opt::operator=(other);
890 rep_opt::operator=(other);
891 matrix_.clear();
892 pivotToColumnIndex_ = other.pivotToColumnIndex_;
893 nextIndex_ = other.nextIndex_;
894 colSettings_ = other.colSettings_;
895
896 matrix_.reserve(other.matrix_.size());
897 for (const auto& cont : other.matrix_){
898 if constexpr (Master_matrix::Option_list::has_map_column_container){
899 _container_insert(cont.second, cont.first);
900 } else {
901 _container_insert(cont);
902 }
903 }
904
905 return *this;
906}
907
908template <class Master_matrix>
909inline void Chain_matrix<Master_matrix>::print() const
910{
911 std::cout << "Column Matrix:\n";
912 if constexpr (!Master_matrix::Option_list::has_map_column_container) {
913 for (id_index i = 0; i < pivotToColumnIndex_.size() && pivotToColumnIndex_[i] != static_cast<index>(-1); ++i) {
914 index pos = pivotToColumnIndex_[i];
915 const Column_type& col = matrix_[pos];
916 for (const auto& cell : col) {
917 std::cout << cell.get_row_index() << " ";
918 }
919 std::cout << "(" << i << ", " << pos << ")\n";
920 }
921 if constexpr (Master_matrix::Option_list::has_row_access) {
922 std::cout << "\n";
923 std::cout << "Row Matrix:\n";
924 for (id_index i = 0; i < pivotToColumnIndex_.size() && pivotToColumnIndex_[i] != static_cast<index>(-1); ++i) {
925 index pos = pivotToColumnIndex_[i];
926 const Row_type& row = ra_opt::get_row(pos);
927 for (const auto& cell : row) {
928 std::cout << cell.get_column_index() << " ";
929 }
930 std::cout << "(" << i << ", " << pos << ")\n";
931 }
932 }
933 } else {
934 for (const auto& p : pivotToColumnIndex_) {
935 const Column_type& col = matrix_.at(p.second);
936 for (const auto& cell : col) {
937 std::cout << cell.get_row_index() << " ";
938 }
939 std::cout << "(" << p.first << ", " << p.second << ")\n";
940 }
941 if constexpr (Master_matrix::Option_list::has_row_access) {
942 std::cout << "\n";
943 std::cout << "Row Matrix:\n";
944 for (const auto& p : pivotToColumnIndex_) {
945 const Row_type& row = ra_opt::get_row(p.first);
946 for (const auto& cell : row) {
947 std::cout << cell.get_column_index() << " ";
948 }
949 std::cout << "(" << p.first << ", " << p.second << ")\n";
950 }
951 }
952 }
953 std::cout << "\n";
954}
955
956template <class Master_matrix>
957template <class Boundary_type>
958inline std::vector<typename Master_matrix::cell_rep_type> Chain_matrix<Master_matrix>::_reduce_boundary(
959 id_index faceID, const Boundary_type& boundary, dimension_type dim)
960{
961 tmp_column_type column(boundary.begin(), boundary.end());
962 if (dim == static_cast<dimension_type>(-1)) dim = boundary.begin() == boundary.end() ? 0 : boundary.size() - 1;
963 std::vector<cell_rep_type> chainsInH; // for corresponding indices in H (paired columns)
964 std::vector<cell_rep_type> chainsInF; // for corresponding indices in F (unpaired, essential columns)
965
966 auto get_last = [&column]() {
967 if constexpr (Master_matrix::Option_list::is_z2)
968 return *(column.rbegin());
969 else
970 return column.rbegin()->first;
971 };
972
973 if (boundary.begin() == boundary.end()) {
974 if constexpr (Master_matrix::Option_list::is_z2)
975 column.insert(faceID);
976 else
977 column.emplace(faceID, 1);
978 _insert_chain(column, dim);
979 return chainsInF;
980 }
981
982 index currentIndex = get_column_with_pivot(get_last());
983
984 while (get_column(currentIndex).is_paired()) {
985 _reduce_by_G(column, chainsInH, currentIndex);
986
987 if (column.empty()) {
988 // produce the sum of all col_h in chains_in_H
989 _build_from_H(faceID, column, chainsInH);
990 // create a new cycle (in F) sigma - \sum col_h
991 _insert_chain(column, dim);
992 return chainsInF;
993 }
994
995 currentIndex = get_column_with_pivot(get_last());
996 }
997
998 while (!column.empty()) {
999 currentIndex = get_column_with_pivot(get_last());
1000
1001 if (!get_column(currentIndex).is_paired()) {
1002 // only fills currentEssentialCycleIndices if Z2 coefficients, so chainsInF remains empty
1003 _reduce_by_F(column, chainsInF, currentIndex);
1004 } else {
1005 _reduce_by_G(column, chainsInH, currentIndex);
1006 }
1007 }
1008
1009 _update_largest_death_in_F(chainsInF);
1010
1011 // Compute the new column zzsh + \sum col_h, for col_h in chains_in_H
1012 _build_from_H(faceID, column, chainsInH);
1013
1014 // Create and insert (\sum col_h) + sigma (in H, paired with chain_fp) in matrix_
1015 if constexpr (Master_matrix::Option_list::is_z2)
1016 _insert_chain(column, dim, chainsInF[0]);
1017 else
1018 _insert_chain(column, dim, chainsInF[0].first);
1019
1020 return chainsInF;
1021}
1022
1023template <class Master_matrix>
1024inline void Chain_matrix<Master_matrix>::_reduce_by_G(tmp_column_type& column,
1025 std::vector<cell_rep_type>& chainsInH,
1026 index currentIndex)
1027{
1028 Column_type& col = get_column(currentIndex);
1029 if constexpr (Master_matrix::Option_list::is_z2) {
1030 _add_to(col, column, 1u); // Reduce with the column col_g
1031 chainsInH.push_back(col.get_paired_chain_index()); // keep the col_h with which col_g is paired
1032 } else {
1033 Field_element_type coef = col.get_pivot_value();
1034 auto& operators = colSettings_->operators;
1035 coef = operators.get_inverse(coef);
1036 operators.multiply_inplace(coef, operators.get_characteristic() - column.rbegin()->second);
1037
1038 _add_to(col, column, coef); // Reduce with the column col_g
1039 chainsInH.emplace_back(col.get_paired_chain_index(), coef); // keep the col_h with which col_g is paired
1040 }
1041}
1042
1043template <class Master_matrix>
1044inline void Chain_matrix<Master_matrix>::_reduce_by_F(tmp_column_type& column,
1045 std::vector<cell_rep_type>& chainsInF,
1046 index currentIndex)
1047{
1048 Column_type& col = get_column(currentIndex);
1049 if constexpr (Master_matrix::Option_list::is_z2) {
1050 _add_to(col, column, 1u); // Reduce with the column col_g
1051 chainsInF.push_back(currentIndex);
1052 } else {
1053 Field_element_type coef = col.get_pivot_value();
1054 auto& operators = colSettings_->operators;
1055 coef = operators.get_inverse(coef);
1056 operators.multiply_inplace(coef, operators.get_characteristic() - column.rbegin()->second);
1057
1058 _add_to(col, column, coef); // Reduce with the column col_g
1059 chainsInF.emplace_back(currentIndex, operators.get_characteristic() - coef);
1060 }
1061}
1062
1063template <class Master_matrix>
1064inline void Chain_matrix<Master_matrix>::_build_from_H(id_index faceID,
1065 tmp_column_type& column,
1066 std::vector<cell_rep_type>& chainsInH)
1067{
1068 if constexpr (Master_matrix::Option_list::is_z2) {
1069 column.insert(faceID);
1070 for (index idx_h : chainsInH) {
1071 _add_to(get_column(idx_h), column, 1u);
1072 }
1073 } else {
1074 column.emplace(faceID, 1);
1075 for (std::pair<index, Field_element_type>& idx_h : chainsInH) {
1076 _add_to(get_column(idx_h.first), column, idx_h.second);
1077 }
1078 }
1079}
1080
1081template <class Master_matrix>
1082inline void Chain_matrix<Master_matrix>::_update_largest_death_in_F(const std::vector<cell_rep_type>& chainsInF)
1083{
1084 if constexpr (Master_matrix::Option_list::is_z2) {
1085 index toUpdate = chainsInF[0];
1086 for (auto other_col_it = chainsInF.begin() + 1; other_col_it != chainsInF.end(); ++other_col_it) {
1087 add_to(*other_col_it, toUpdate);
1088 }
1089 } else {
1090 index toUpdate = chainsInF[0].first;
1091 get_column(toUpdate) *= chainsInF[0].second;
1092 for (auto other_col_it = chainsInF.begin() + 1; other_col_it != chainsInF.end(); ++other_col_it) {
1093 multiply_source_and_add_to(other_col_it->second, other_col_it->first, toUpdate);
1094 }
1095 }
1096}
1097
1098template <class Master_matrix>
1099inline void Chain_matrix<Master_matrix>::_insert_chain(const tmp_column_type& column, dimension_type dimension)
1100{
1101 _container_insert(column, nextIndex_, dimension);
1102 _add_bar(dimension);
1103
1104 ++nextIndex_;
1105}
1106
1107template <class Master_matrix>
1108inline void Chain_matrix<Master_matrix>::_insert_chain(const tmp_column_type& column,
1109 dimension_type dimension,
1110 index pair)
1111{
1112 // true when no vine updates and if nextIndex_ is updated in remove_last for special case of no vines
1113 // because then @ref PosIdx == @ref MatIdx
1114 pos_index pairPos = pair;
1115
1116 _container_insert(column, nextIndex_, dimension);
1117
1118 get_column(nextIndex_).assign_paired_chain(pair);
1119 auto& pairCol = get_column(pair);
1120 pairCol.assign_paired_chain(nextIndex_);
1121
1122 if constexpr (Master_matrix::Option_list::has_column_pairings && Master_matrix::Option_list::has_vine_update) {
1123 pairPos = swap_opt::CP::pivotToPosition_[pairCol.get_pivot()];
1124 }
1125
1126 _update_barcode(pairPos);
1127
1128 ++nextIndex_;
1129}
1130
1131template <class Master_matrix>
1132inline void Chain_matrix<Master_matrix>::_add_to(const Column_type& column,
1133 tmp_column_type& set,
1134 [[maybe_unused]] unsigned int coef)
1135{
1136 if constexpr (Master_matrix::Option_list::is_z2) {
1137 std::pair<typename std::set<index>::iterator, bool> res_insert;
1138 for (const Cell& cell : column) {
1139 res_insert = set.insert(cell.get_row_index());
1140 if (!res_insert.second) {
1141 set.erase(res_insert.first);
1142 }
1143 }
1144 } else {
1145 auto& operators = colSettings_->operators;
1146 for (const Cell& cell : column) {
1147 auto res = set.emplace(cell.get_row_index(), cell.get_element());
1148 if (res.second){
1149 operators.multiply_inplace(res.first->second, coef);
1150 } else {
1151 operators.multiply_and_add_inplace_back(cell.get_element(), coef, res.first->second);
1152 if (res.first->second == Field_operators::get_additive_identity()) {
1153 set.erase(res.first);
1154 }
1155 }
1156 }
1157 }
1158}
1159
1160template <class Master_matrix>
1161template <typename F>
1162inline void Chain_matrix<Master_matrix>::_add_to(Column_type& target, F&& addition)
1163{
1164 auto pivot = target.get_pivot();
1165 addition();
1166
1167 if (pivot != target.get_pivot()) {
1168 if constexpr (Master_matrix::Option_list::has_map_column_container) {
1169 std::swap(pivotToColumnIndex_.at(pivot), pivotToColumnIndex_.at(target.get_pivot()));
1170 } else {
1171 std::swap(pivotToColumnIndex_[pivot], pivotToColumnIndex_[target.get_pivot()]);
1172 }
1173 }
1174}
1175
1176template <class Master_matrix>
1177inline void Chain_matrix<Master_matrix>::_remove_last(index lastIndex)
1178{
1179 static_assert(Master_matrix::Option_list::has_removable_columns,
1180 "'_remove_last' is not implemented for the chosen options.");
1181 static_assert(Master_matrix::Option_list::has_map_column_container || !Master_matrix::Option_list::has_vine_update,
1182 "'_remove_last' is not implemented for the chosen options.");
1183
1184 id_index pivot;
1185
1186 if constexpr (Master_matrix::Option_list::has_map_column_container) {
1187 auto itToErase = matrix_.find(lastIndex);
1188 Column_type& colToErase = itToErase->second;
1189 pivot = colToErase.get_pivot();
1190
1191 if constexpr (Master_matrix::Option_list::has_matrix_maximal_dimension_access) {
1192 dim_opt::update_down(colToErase.get_dimension());
1193 }
1194
1195 if (colToErase.is_paired()) matrix_.at(colToErase.get_paired_chain_index()).unassign_paired_chain();
1196 pivotToColumnIndex_.erase(pivot);
1197 matrix_.erase(itToErase);
1198 } else {
1199 GUDHI_CHECK(lastIndex == nextIndex_ - 1 && nextIndex_ == matrix_.size(),
1200 std::logic_error("Chain_matrix::_remove_last - Indexation problem."));
1201
1202 Column_type& colToErase = matrix_[lastIndex];
1203 pivot = colToErase.get_pivot();
1204
1205 if constexpr (Master_matrix::Option_list::has_matrix_maximal_dimension_access) {
1206 dim_opt::update_down(colToErase.get_dimension());
1207 }
1208
1209 if (colToErase.is_paired()) matrix_.at(colToErase.get_paired_chain_index()).unassign_paired_chain();
1210 pivotToColumnIndex_[pivot] = -1;
1211 matrix_.pop_back();
1212 // TODO: resize matrix_ when a lot is removed? Could be not the best strategy if user inserts a lot back afterwards.
1213 }
1214
1215 if constexpr (!Master_matrix::Option_list::has_vine_update) {
1216 --nextIndex_; // should not be updated when there are vine updates, as possibly lastIndex != nextIndex - 1
1217 }
1218
1219 if constexpr (Master_matrix::Option_list::has_column_pairings) {
1220 auto it = _indexToBar().find(--_nextPosition());
1221 typename barcode_type::iterator bar = it->second;
1222
1223 if (bar->death == static_cast<pos_index>(-1))
1224 _barcode().erase(bar);
1225 else
1226 bar->death = -1;
1227
1228 _indexToBar().erase(it);
1229 if constexpr (Master_matrix::Option_list::has_vine_update) swap_opt::CP::pivotToPosition_.erase(pivot);
1230 }
1231
1232 if constexpr (Master_matrix::Option_list::has_row_access) {
1233 GUDHI_CHECK(
1234 ra_opt::get_row(pivot).size() == 0,
1235 std::invalid_argument(
1236 "Chain_matrix::_remove_last - Column asked to be removed does not corresponds to a maximal simplex."));
1237 if constexpr (Master_matrix::Option_list::has_removable_rows) {
1238 ra_opt::erase_empty_row(pivot);
1239 }
1240 }
1241}
1242
1243template <class Master_matrix>
1244inline void Chain_matrix<Master_matrix>::_update_barcode(pos_index birth)
1245{
1246 if constexpr (Master_matrix::Option_list::has_column_pairings) {
1247 if constexpr (Master_matrix::Option_list::has_removable_columns) {
1248 auto& barIt = _indexToBar().at(birth);
1249 barIt->death = _nextPosition();
1250 _indexToBar().try_emplace(_nextPosition(), barIt); // list so iterators are stable
1251 } else {
1252 _barcode()[_indexToBar()[birth]].death = _nextPosition();
1253 _indexToBar().push_back(_indexToBar()[birth]);
1254 }
1255 ++_nextPosition();
1256 }
1257}
1258
1259template <class Master_matrix>
1260inline void Chain_matrix<Master_matrix>::_add_bar(dimension_type dim)
1261{
1262 if constexpr (Master_matrix::Option_list::has_column_pairings) {
1263 _barcode().emplace_back(dim, _nextPosition(), -1);
1264 if constexpr (Master_matrix::Option_list::has_removable_columns) {
1265 _indexToBar().try_emplace(_nextPosition(), --_barcode().end());
1266 } else {
1267 _indexToBar().push_back(_barcode().size() - 1);
1268 }
1269 ++_nextPosition();
1270 }
1271}
1272
1273template <class Master_matrix>
1274template <class Container_type>
1275inline void Chain_matrix<Master_matrix>::_container_insert(const Container_type& column,
1276 index pos,
1277 dimension_type dim)
1278{
1279 id_index pivot;
1280 if constexpr (Master_matrix::Option_list::is_z2) {
1281 pivot = *(column.rbegin());
1282 } else {
1283 pivot = column.rbegin()->first;
1284 }
1285 if constexpr (Master_matrix::Option_list::has_map_column_container) {
1286 pivotToColumnIndex_.try_emplace(pivot, pos);
1287 if constexpr (Master_matrix::Option_list::has_row_access) {
1288 matrix_.try_emplace(pos, Column_type(pos, column, dim, ra_opt::rows_, colSettings_));
1289 } else {
1290 matrix_.try_emplace(pos, Column_type(column, dim, colSettings_));
1291 }
1292 } else {
1293 if constexpr (Master_matrix::Option_list::has_row_access) {
1294 matrix_.emplace_back(pos, column, dim, ra_opt::rows_, colSettings_);
1295 } else {
1296 matrix_.emplace_back(column, dim, colSettings_);
1297 }
1298 pivotToColumnIndex_[pivot] = pos;
1299 }
1300}
1301
1302template <class Master_matrix>
1303inline void Chain_matrix<Master_matrix>::_container_insert(const Column_type& column, [[maybe_unused]] index pos)
1304{
1305 if constexpr (Master_matrix::Option_list::has_map_column_container) {
1306 if constexpr (Master_matrix::Option_list::has_row_access) {
1307 matrix_.try_emplace(pos, Column_type(column, column.get_column_index(), ra_opt::rows_, colSettings_));
1308 } else {
1309 matrix_.try_emplace(pos, Column_type(column, colSettings_));
1310 }
1311 } else {
1312 if constexpr (Master_matrix::Option_list::has_row_access) {
1313 matrix_.emplace_back(column, column.get_column_index(), ra_opt::rows_, colSettings_);
1314 } else {
1315 matrix_.emplace_back(column, colSettings_);
1316 }
1317 }
1318}
1319
1320template <class Master_matrix>
1322{
1323 if constexpr (Master_matrix::Option_list::has_vine_update)
1324 return swap_opt::template Chain_pairing<Master_matrix>::barcode_;
1325 else
1326 return pair_opt::barcode_;
1327}
1328
1329template <class Master_matrix>
1330inline constexpr typename Chain_matrix<Master_matrix>::bar_dictionnary_type&
1332{
1333 if constexpr (Master_matrix::Option_list::has_vine_update)
1334 return swap_opt::template Chain_pairing<Master_matrix>::indexToBar_;
1335 else
1336 return pair_opt::indexToBar_;
1337}
1338
1339template <class Master_matrix>
1341{
1342 if constexpr (Master_matrix::Option_list::has_vine_update)
1343 return swap_opt::template Chain_pairing<Master_matrix>::nextPosition_;
1344 else
1345 return pair_opt::nextPosition_;
1346}
1347
1348} // namespace persistence_matrix
1349} // namespace Gudhi
1350
1351#endif // PM_CHAIN_MATRIX_H
Matrix structure storing a compatible base of a filtered chain complex. See . The base is constructed...
Definition chain_matrix.h:50
typename Master_matrix::Cell_constructor Cell_constructor
Definition chain_matrix.h:61
friend void swap(Chain_matrix &matrix1, Chain_matrix &matrix2)
Swap operator.
Definition chain_matrix.h:452
void remove_maximal_face(id_index faceID)
Only available if PersistenceMatrixOptions::has_removable_columns and PersistenceMatrixOptions::has_v...
Definition chain_matrix.h:728
typename Master_matrix::Cell_type Cell
Definition chain_matrix.h:60
std::vector< cell_rep_type > insert_boundary(id_index faceID, const Boundary_type &boundary, dimension_type dim=-1)
It does the same as the other version, but allows the boundary faces to be identified without restric...
void reset(Column_settings *colSettings)
Resets the matrix to an empty matrix.
Definition chain_matrix.h:438
typename Master_matrix::pos_index pos_index
Definition chain_matrix.h:68
id_index get_pivot(index columnIndex)
Returns the row index of the pivot of the given column.
Definition chain_matrix.h:879
dimension_type get_column_dimension(index columnIndex) const
Returns the dimension of the given column.
Definition chain_matrix.h:824
typename Master_matrix::cell_rep_type cell_rep_type
Definition chain_matrix.h:65
typename Master_matrix::boundary_type boundary_type
Definition chain_matrix.h:64
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:862
typename Master_matrix::Row_type Row_type
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::index index
Definition chain_matrix.h:66
index get_number_of_columns() const
Returns the current number of columns in the matrix.
Definition chain_matrix.h:814
typename Master_matrix::element_type Field_element_type
Definition chain_matrix.h:56
Column_type & get_column(index columnIndex)
Returns the column at the given MatIdx index. The type of the column depends on the choosen options,...
Definition chain_matrix.h:707
void multiply_source_and_add_to(const Field_element_type &coefficient, index sourceColumnIndex, index targetColumnIndex)
Multiplies the source column with the coefficiant before adding it to the target column....
Definition chain_matrix.h:847
void multiply_target_and_add_to(index sourceColumnIndex, const Field_element_type &coefficient, index targetColumnIndex)
Multiplies the target column with the coefficiant and then adds the source column to it....
Definition chain_matrix.h:838
typename Master_matrix::dimension_type dimension_type
Definition chain_matrix.h:69
void remove_last()
Only available if PersistenceMatrixOptions::has_removable_columns is true and, if PersistenceMatrixOp...
Definition chain_matrix.h:785
std::vector< cell_rep_type > insert_boundary(const Boundary_type &boundary, dimension_type dim=-1)
Inserts at the end of the matrix a new ordered column corresponding to the given boundary....
index get_column_with_pivot(id_index faceID) const
Returns the column with given row index as pivot. Assumes that the pivot exists.
Definition chain_matrix.h:868
void add_to(index sourceColumnIndex, index targetColumnIndex)
Adds column at sourceColumnIndex onto the column at targetColumnIndex in the matrix.
Definition chain_matrix.h:831
typename Master_matrix::id_index id_index
Definition chain_matrix.h:67
bool is_zero_cell(index columnIndex, id_index rowIndex) const
Indicates if the cell at given coordinates has value zero.
Definition chain_matrix.h:856
Chain_matrix & operator=(const Chain_matrix &other)
Assign operator.
Definition chain_matrix.h:885
typename Master_matrix::Column_settings Column_settings
Definition chain_matrix.h:63
typename Master_matrix::Column_type Column_type
Definition chain_matrix.h:57
Overlay for non-basic matrices replacing all input and output MatIdx indices of the original methods ...
Definition overlay_ididx_to_matidx.h:42
Data structure for matrices, and in particular thought for matrices representing filtered complexes i...
Definition matrix.h:143
bool is_zero_cell(index columnIndex, id_index rowIndex)
Indicates if the cell at given coordinates has value zero.
Definition matrix.h:1880
bool is_zero_column(index columnIndex)
Indicates if the column at given index has value zero.
Definition matrix.h:1899
returned_column_type & get_column(index columnIndex)
Returns the column at the given MatIdx index. For RU matrices, is equivalent to get_column(columnInde...
Definition matrix.h:1609
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
void remove_maximal_face(index columnIndex)
Only available for RU and chain matrices and if PersistenceMatrixOptions::has_removable_columns and P...
Definition matrix.h:1690
typename PersistenceMatrixOptions::index_type pos_index
Definition matrix.h:148
typename std::conditional< hasFixedBarcode, std::vector< Bar >, typename std::conditional< PersistenceMatrixOptions::has_removable_columns, std::list< Bar >, std::vector< Bar > >::type >::type barcode_type
Type of the computed barcode. It is either a list of Matrix::Bar or a vector of Matrix::Bar,...
Definition matrix.h:450
dimension_type get_column_dimension(index columnIndex) const
Returns the dimension of the given face. Only available for non-basic matrices.
Definition matrix.h:1747
typename std::conditional< PersistenceMatrixOptions::column_type==Column_types::HEAP, Heap_column_type, typename std::conditional< PersistenceMatrixOptions::column_type==Column_types::LIST, List_column_type, typename std::conditional< PersistenceMatrixOptions::column_type==Column_types::SET, Set_column_type, typename std::conditional< PersistenceMatrixOptions::column_type==Column_types::UNORDERED_SET, Unordered_set_column_type, typename std::conditional< PersistenceMatrixOptions::column_type==Column_types::VECTOR, Vector_column_type, typename std::conditional< PersistenceMatrixOptions::column_type==Column_types::INTRUSIVE_LIST, Intrusive_list_column_type, typename std::conditional< PersistenceMatrixOptions::column_type==Column_types::NAIVE_VECTOR, Naive_vector_column_type, Intrusive_set_column_type >::type >::type >::type >::type >::type >::type >::type Column_type
Type of the columns stored in the matrix. The type depends on the value of PersistenceMatrixOptions::...
Definition matrix.h:370
typename PersistenceMatrixOptions::index_type index
Definition matrix.h:146
std::enable_if_t< std::is_integral_v< Index_type > > add_to(Index_type sourceColumnIndex, Index_type targetColumnIndex)
Adds column at sourceColumnIndex onto the column at targetColumnIndex in the matrix....
Definition matrix.h:1757
insertion_return_type insert_boundary(const Boundary_type &boundary, dimension_type dim=-1)
Inserts at the end of the matrix a new ordered column corresponding to the given boundary....
Definition matrix.h:1572
std::enable_if_t< std::is_integral_v< Index_type > > multiply_target_and_add_to(Index_type sourceColumnIndex, int coefficient, Index_type targetColumnIndex)
Multiplies the target column with the coefficient and then adds the source column to it....
Definition matrix.h:1777
void remove_last()
Removes the last inserted column/face from the matrix. If the matrix is non basic,...
Definition matrix.h:1718
index get_column_with_pivot(id_index faceIndex) const
Returns the MatIdx index of the column which has the given row index as pivot. Only available for RU ...
Definition matrix.h:1918
typename PersistenceMatrixOptions::index_type id_index
Definition matrix.h:147
Matrix & operator=(Matrix other)
Assign operator.
Definition matrix.h:1939
std::enable_if_t< std::is_integral_v< Index_type > > multiply_source_and_add_to(int coefficient, Index_type sourceColumnIndex, Index_type targetColumnIndex)
Multiplies the source column with the coefficient before adding it to the target column....
Definition matrix.h:1807
index get_number_of_columns() const
Returns the current number of columns in the matrix.
Definition matrix.h:1741
typename PersistenceMatrixOptions::dimension_type dimension_type
Definition matrix.h:149
Gudhi namespace.
Definition SimplicialComplexForAlpha.h:14
Contains the Id_to_index_overlay class.