Loading...
Searching...
No Matches
column_utilities.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) 2024 Inria
6 *
7 * Modification(s):
8 * - YYYY/MM Author: Description of the modification
9 */
10
17#ifndef PM_COLUMN_UTILITIES_H
18#define PM_COLUMN_UTILITIES_H
19
20#include <cstddef>
21#include <stdexcept>
22
24
25namespace Gudhi {
26namespace persistence_matrix {
27
28template <class Cell, typename Cell_iterator>
29Cell* _get_cell(const Cell_iterator& itTarget)
30{
31 if constexpr (Cell::Master::Option_list::column_type == Column_types::INTRUSIVE_LIST ||
32 Cell::Master::Option_list::column_type == Column_types::INTRUSIVE_SET) {
33 return &*itTarget;
34 } else {
35 return *itTarget;
36 }
37}
38
39// works only for ordered columns
40template <class Column_type, class Cell_iterator, typename F1, typename F2, typename F3, typename F4>
41void _generic_merge_cell_to_column(Column_type& targetColumn,
42 Cell_iterator& itSource,
43 typename Column_type::Column_type::iterator& itTarget,
44 F1&& process_target,
45 F2&& process_source,
46 F3&& update_target1,
47 F4&& update_target2,
48 bool& pivotIsZeroed)
49{
50 typename Column_type::Cell* targetCell = _get_cell<typename Column_type::Cell>(itTarget);
51
52 if (targetCell->get_row_index() < itSource->get_row_index()) {
53 process_target(targetCell);
54 ++itTarget;
55 } else if (targetCell->get_row_index() > itSource->get_row_index()) {
56 process_source(itSource, itTarget);
57 ++itSource;
58 } else {
59 if constexpr (Column_type::Master::Option_list::is_z2) {
60 //_multiply_*_and_add never enters here so not treated
61 if constexpr (Column_type::Master::isNonBasic && !Column_type::Master::Option_list::is_of_boundary_type) {
62 if (targetCell->get_row_index() == targetColumn.get_pivot()) pivotIsZeroed = true;
63 }
64 targetColumn._delete_cell(itTarget);
65 } else {
66 update_target1(targetCell->get_element(), itSource);
67 if (targetCell->get_element() == Column_type::Field_operators::get_additive_identity()) {
68 if constexpr (Column_type::Master::isNonBasic && !Column_type::Master::Option_list::is_of_boundary_type) {
69 if (targetCell->get_row_index() == targetColumn.get_pivot()) pivotIsZeroed = true;
70 }
71 targetColumn._delete_cell(itTarget);
72 } else {
73 update_target2(targetCell);
74 if constexpr (Column_type::Master::Option_list::has_row_access)
75 targetColumn.update_cell(*targetCell);
76 ++itTarget;
77 }
78 }
79 ++itSource;
80 }
81}
82
83// works only for ordered columns
84template <class Column_type, class Cell_range, typename F1, typename F2, typename F3, typename F4, typename F5>
85bool _generic_add_to_column(const Cell_range& source,
86 Column_type& targetColumn,
87 F1&& process_target,
88 F2&& process_source,
89 F3&& update_target1,
90 F4&& update_target2,
91 F5&& finish_target)
92{
93 bool pivotIsZeroed = false;
94
95 auto& target = targetColumn.column_;
96 auto itTarget = target.begin();
97 auto itSource = source.begin();
98 while (itTarget != target.end() && itSource != source.end()) {
99 _generic_merge_cell_to_column(targetColumn, itSource, itTarget,
100 process_target, process_source, update_target1, update_target2,
101 pivotIsZeroed);
102 }
103
104 finish_target(itTarget);
105
106 while (itSource != source.end()) {
107 process_source(itSource, target.end());
108 ++itSource;
109 }
110
111 return pivotIsZeroed;
112}
113
114template <class Column_type, class Cell_range>
115bool _add_to_column(const Cell_range& source, Column_type& targetColumn)
116{
117 return _generic_add_to_column(
118 source,
119 targetColumn, [&]([[maybe_unused]] typename Column_type::Cell* cellTarget) {},
120 [&](typename Cell_range::const_iterator& itSource, const typename Column_type::Column_type::iterator& itTarget) {
121 if constexpr (Column_type::Master::Option_list::is_z2) {
122 targetColumn._insert_cell(itSource->get_row_index(), itTarget);
123 } else {
124 targetColumn._insert_cell(itSource->get_element(), itSource->get_row_index(), itTarget);
125 }
126 },
127 [&](typename Column_type::Field_element_type& targetElement, typename Cell_range::const_iterator& itSource) {
128 if constexpr (!Column_type::Master::Option_list::is_z2)
129 targetColumn.operators_->add_inplace(targetElement, itSource->get_element());
130 },
131 [&]([[maybe_unused]] typename Column_type::Cell* cellTarget) {},
132 [&]([[maybe_unused]] typename Column_type::Column_type::iterator& itTarget) {});
133}
134
135template <class Column_type, class Cell_range>
136bool _multiply_target_and_add_to_column(const typename Column_type::Field_element_type& val,
137 const Cell_range& source,
138 Column_type& targetColumn)
139{
140 if (val == 0u) {
141 if constexpr (Column_type::Master::isNonBasic && !Column_type::Master::Option_list::is_of_boundary_type) {
142 throw std::invalid_argument("A chain column should not be multiplied by 0.");
143 // this would not only mess up the base, but also the pivots stored.
144 } else {
145 targetColumn.clear();
146 }
147 }
148
149 return _generic_add_to_column(
150 source,
151 targetColumn,
152 [&](typename Column_type::Cell* cellTarget) {
153 targetColumn.operators_->multiply_inplace(cellTarget->get_element(), val);
154 // targetColumn.ra_opt::update_cell(*itTarget) produces an internal compiler error
155 // even though it works in _generic_add_to_column... Probably because of the lambda.
156 if constexpr (Column_type::Master::Option_list::has_row_access)
157 targetColumn.update_cell(*cellTarget);
158 },
159 [&](typename Cell_range::const_iterator& itSource, const typename Column_type::Column_type::iterator& itTarget) {
160 targetColumn._insert_cell(itSource->get_element(), itSource->get_row_index(), itTarget);
161 },
162 [&](typename Column_type::Field_element_type& targetElement, typename Cell_range::const_iterator& itSource) {
163 targetColumn.operators_->multiply_and_add_inplace_front(targetElement, val, itSource->get_element());
164 },
165 [&]([[maybe_unused]] typename Column_type::Cell* cellTarget) {},
166 [&](typename Column_type::Column_type::iterator& itTarget) {
167 while (itTarget != targetColumn.column_.end()) {
168 typename Column_type::Cell* targetCell = _get_cell<typename Column_type::Cell>(itTarget);
169 targetColumn.operators_->multiply_inplace(targetCell->get_element(), val);
170 if constexpr (Column_type::Master::Option_list::has_row_access)
171 targetColumn.update_cell(*targetCell);
172 itTarget++;
173 }
174 });
175}
176
177template <class Column_type, class Cell_range>
178bool _multiply_source_and_add_to_column(const typename Column_type::Field_element_type& val, const Cell_range& source,
179 Column_type& targetColumn)
180{
181 if (val == 0u) {
182 return false;
183 }
184
185 return _generic_add_to_column(
186 source,
187 targetColumn, []([[maybe_unused]] typename Column_type::Cell* cellTarget) {},
188 [&](typename Cell_range::const_iterator& itSource, const typename Column_type::Column_type::iterator& itTarget) {
189 typename Column_type::Cell* cell =
190 targetColumn._insert_cell(itSource->get_element(), itSource->get_row_index(), itTarget);
191 targetColumn.operators_->multiply_inplace(cell->get_element(), val);
192 if constexpr (Column_type::Master::Option_list::has_row_access)
193 targetColumn.update_cell(*cell);
194 },
195 [&](typename Column_type::Field_element_type& targetElement, typename Cell_range::const_iterator& itSource) {
196 targetColumn.operators_->multiply_and_add_inplace_back(itSource->get_element(), val, targetElement);
197 },
198 [&]([[maybe_unused]] typename Column_type::Cell* cellTarget) {},
199 []([[maybe_unused]] typename Column_type::Column_type::iterator& itTarget) {});
200}
201
202// column has to be ordered (ie. not suited for unordered_map and heap) and contain the exact values
203// (ie. not suited for vector and heap). A same colonne but ordered differently will have another hash value.
204template <class Column_type>
205std::size_t hash_column(const Column_type& column) {
206 std::size_t seed = 0;
207 for (auto& cell : column) {
208 seed ^= std::hash<unsigned int>()(cell.get_row_index() * static_cast<unsigned int>(cell.get_element())) +
209 0x9e3779b9 + (seed << 6) + (seed >> 2);
210 }
211 return seed;
212}
213
214} // namespace persistence_matrix
215} // namespace Gudhi
216
217#endif // PM_COLUMN_UTILITIES_H
@ INTRUSIVE_LIST
Definition persistence_matrix_options.h:39
@ INTRUSIVE_SET
Definition persistence_matrix_options.h:40
Gudhi namespace.
Definition SimplicialComplexForAlpha.h:14
Contains the options for the matrix template.