Loading...
Searching...
No Matches
vector_column.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
18#ifndef PM_VECTOR_COLUMN_H
19#define PM_VECTOR_COLUMN_H
20
21#include <cstddef>
22#include <vector>
23#include <stdexcept>
24#include <type_traits>
25#include <algorithm> //binary_search
26#include <unordered_set>
27#include <utility> //std::swap, std::move & std::exchange
28
29#include <boost/iterator/indirect_iterator.hpp>
31
32namespace Gudhi {
33namespace persistence_matrix {
34
49template <class Master_matrix>
50class Vector_column : public Master_matrix::Row_access_option,
51 public Master_matrix::Column_dimension_option,
52 public Master_matrix::Chain_column_option
53{
54 public:
55 using Master = Master_matrix;
56 using index = typename Master_matrix::index;
57 using id_index = typename Master_matrix::id_index;
58 using dimension_type = typename Master_matrix::dimension_type;
59 using Field_element_type = typename Master_matrix::element_type;
60 using Cell = typename Master_matrix::Cell_type;
61 using Column_settings = typename Master_matrix::Column_settings;
62
63 private:
64 using Field_operators = typename Master_matrix::Field_operators;
65 using Column_type = std::vector<Cell*>;
66 using Cell_constructor = typename Master_matrix::Cell_constructor;
67
68 public:
69 using iterator = boost::indirect_iterator<typename Column_type::iterator>;
70 using const_iterator = boost::indirect_iterator<typename Column_type::const_iterator>;
71 using reverse_iterator = boost::indirect_iterator<typename Column_type::reverse_iterator>;
72 using const_reverse_iterator = boost::indirect_iterator<typename Column_type::const_reverse_iterator>;
73
74 Vector_column(Column_settings* colSettings = nullptr);
75 template <class Container_type = typename Master_matrix::boundary_type>
76 Vector_column(const Container_type& nonZeroRowIndices, Column_settings* colSettings);
77 template <class Container_type = typename Master_matrix::boundary_type, class Row_container_type>
78 Vector_column(index columnIndex,
79 const Container_type& nonZeroRowIndices,
80 Row_container_type* rowContainer,
81 Column_settings* colSettings);
82 template <class Container_type = typename Master_matrix::boundary_type>
83 Vector_column(const Container_type& nonZeroChainRowIndices,
84 dimension_type dimension,
85 Column_settings* colSettings);
86 template <class Container_type = typename Master_matrix::boundary_type, class Row_container_type>
87 Vector_column(index columnIndex,
88 const Container_type& nonZeroChainRowIndices,
89 dimension_type dimension,
90 Row_container_type* rowContainer,
91 Column_settings* colSettings);
92 Vector_column(const Vector_column& column,
93 Column_settings* colSettings = nullptr);
94 template <class Row_container_type>
95 Vector_column(const Vector_column& column,
96 index columnIndex,
97 Row_container_type* rowContainer,
98 Column_settings* colSettings = nullptr);
99 Vector_column(Vector_column&& column) noexcept;
101
102 std::vector<Field_element_type> get_content(int columnLength = -1) const;
103 bool is_non_zero(id_index rowIndex) const;
104 bool is_empty() const;
105 std::size_t size() const;
106
107 template <class Map_type>
108 void reorder(const Map_type& valueMap, [[maybe_unused]] index columnIndex = -1);
109 void clear();
110 // do not clear a cell to 0 if the cell was already 0, otherwise size/is_empty will be wrong.
111 void clear(id_index rowIndex);
112
113 id_index get_pivot();
114 Field_element_type get_pivot_value();
115
116 iterator begin() noexcept;
117 const_iterator begin() const noexcept;
118 iterator end() noexcept;
119 const_iterator end() const noexcept;
120 reverse_iterator rbegin() noexcept;
121 const_reverse_iterator rbegin() const noexcept;
122 reverse_iterator rend() noexcept;
123 const_reverse_iterator rend() const noexcept;
124
125 template <class Cell_range>
126 Vector_column& operator+=(const Cell_range& column);
127 Vector_column& operator+=(Vector_column& column);
128
129 Vector_column& operator*=(unsigned int v);
130
131 // this = v * this + column
132 template <class Cell_range>
133 Vector_column& multiply_target_and_add(const Field_element_type& val, const Cell_range& column);
134 Vector_column& multiply_target_and_add(const Field_element_type& val, Vector_column& column);
135 // this = this + column * v
136 template <class Cell_range>
137 Vector_column& multiply_source_and_add(const Cell_range& column, const Field_element_type& val);
138 Vector_column& multiply_source_and_add(Vector_column& column, const Field_element_type& val);
139
140 std::size_t compute_hash_value();
141
142 friend bool operator==(const Vector_column& c1, const Vector_column& c2) {
143 if (&c1 == &c2) return true;
144 if (c1.erasedValues_.empty() && c2.erasedValues_.empty() && c1.column_.size() != c2.column_.size()) return false;
145
146 auto it1 = c1.column_.begin();
147 auto it2 = c2.column_.begin();
148 while (it1 != c1.column_.end() && it2 != c2.column_.end()) {
149 if constexpr (!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type) {
150 while (it1 != c1.column_.end() && c1.erasedValues_.find((*it1)->get_row_index()) != c1.erasedValues_.end())
151 ++it1;
152 while (it2 != c2.column_.end() && c2.erasedValues_.find((*it2)->get_row_index()) != c2.erasedValues_.end())
153 ++it2;
154 if (it1 == c1.column_.end() || it2 == c2.column_.end()) break;
155 }
156 if constexpr (Master_matrix::Option_list::is_z2) {
157 if ((*it1)->get_row_index() != (*it2)->get_row_index()) return false;
158 } else {
159 if ((*it1)->get_row_index() != (*it2)->get_row_index() || (*it1)->get_element() != (*it2)->get_element())
160 return false;
161 }
162 ++it1;
163 ++it2;
164 }
165
166 if constexpr (!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type) {
167 while (it1 != c1.column_.end() && c1.erasedValues_.find((*it1)->get_row_index()) != c1.erasedValues_.end()) ++it1;
168 while (it2 != c2.column_.end() && c2.erasedValues_.find((*it2)->get_row_index()) != c2.erasedValues_.end()) ++it2;
169 return it2 == c2.column_.end() && it1 == c1.column_.end();
170 } else {
171 return true;
172 }
173 }
174 friend bool operator<(const Vector_column& c1, const Vector_column& c2) {
175 if (&c1 == &c2) return false;
176
177 auto it1 = c1.column_.begin();
178 auto it2 = c2.column_.begin();
179 while (it1 != c1.column_.end() && it2 != c2.column_.end()) {
180 if constexpr (!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type) {
181 while (it1 != c1.column_.end() && c1.erasedValues_.find((*it1)->get_row_index()) != c1.erasedValues_.end())
182 ++it1;
183 while (it2 != c2.column_.end() && c2.erasedValues_.find((*it2)->get_row_index()) != c2.erasedValues_.end())
184 ++it2;
185 if (it1 == c1.column_.end() || it2 == c2.column_.end()) break;
186 }
187
188 if ((*it1)->get_row_index() != (*it2)->get_row_index()) return (*it1)->get_row_index() < (*it2)->get_row_index();
189 if constexpr (!Master_matrix::Option_list::is_z2) {
190 if ((*it1)->get_element() != (*it2)->get_element()) return (*it1)->get_element() < (*it2)->get_element();
191 }
192 ++it1;
193 ++it2;
194 }
195 if constexpr (!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type) {
196 while (it1 != c1.column_.end() && c1.erasedValues_.find((*it1)->get_row_index()) != c1.erasedValues_.end()) ++it1;
197 while (it2 != c2.column_.end() && c2.erasedValues_.find((*it2)->get_row_index()) != c2.erasedValues_.end()) ++it2;
198 }
199 return it2 != c2.column_.end();
200 }
201
202 // Disabled with row access.
203 Vector_column& operator=(const Vector_column& other);
204
205 friend void swap(Vector_column& col1, Vector_column& col2) {
206 swap(static_cast<typename Master_matrix::Row_access_option&>(col1),
207 static_cast<typename Master_matrix::Row_access_option&>(col2));
208 swap(static_cast<typename Master_matrix::Column_dimension_option&>(col1),
209 static_cast<typename Master_matrix::Column_dimension_option&>(col2));
210 swap(static_cast<typename Master_matrix::Chain_column_option&>(col1),
211 static_cast<typename Master_matrix::Chain_column_option&>(col2));
212 col1.column_.swap(col2.column_);
213 col1.erasedValues_.swap(col2.erasedValues_);
214 std::swap(col1.operators_, col2.operators_);
215 std::swap(col1.cellPool_, col2.cellPool_);
216 }
217
218 private:
219 using ra_opt = typename Master_matrix::Row_access_option;
220 using dim_opt = typename Master_matrix::Column_dimension_option;
221 using chain_opt = typename Master_matrix::Chain_column_option;
222
223 Column_type column_;
224 std::unordered_set<id_index> erasedValues_; // TODO: test other containers? Useless when clear(index) is never
225 // called, how much is it worth it?
226 Field_operators* operators_;
227 Cell_constructor* cellPool_;
228
229 template <class Column_type, class Cell_iterator, typename F1, typename F2, typename F3, typename F4>
230 friend void _generic_merge_cell_to_column(Column_type& targetColumn,
231 Cell_iterator& itSource,
232 typename Column_type::Column_type::iterator& itTarget,
233 F1&& process_target,
234 F2&& process_source,
235 F3&& update_target1,
236 F4&& update_target2,
237 bool& pivotIsZeroed);
238
239 void _delete_cell(Cell* cell);
240 void _delete_cell(typename Column_type::iterator& it);
241 Cell* _insert_cell(const Field_element_type& value, id_index rowIndex, Column_type& column);
242 void _insert_cell(id_index rowIndex, Column_type& column);
243 void _update_cell(const Field_element_type& value, id_index rowIndex, index position);
244 void _update_cell(id_index rowIndex, index position);
245 template <class Cell_range>
246 bool _add(const Cell_range& column);
247 template <class Cell_range>
248 bool _multiply_target_and_add(const Field_element_type& val, const Cell_range& column);
249 template <class Cell_range>
250 bool _multiply_source_and_add(const Cell_range& column, const Field_element_type& val);
251 template <class Cell_range, typename F1, typename F2, typename F3, typename F4>
252 bool _generic_add(const Cell_range& source,
253 F1&& process_target,
254 F2&& process_source,
255 F3&& update_target1,
256 F4&& update_target2);
257};
258
259template <class Master_matrix>
260inline Vector_column<Master_matrix>::Vector_column(Column_settings* colSettings)
261 : ra_opt(), dim_opt(), chain_opt(), operators_(nullptr), cellPool_(colSettings == nullptr ? nullptr : &(colSettings->cellConstructor))
262{
263 if (operators_ == nullptr && cellPool_ == nullptr) return; //to allow default constructor which gives a dummy column
264 if constexpr (!Master_matrix::Option_list::is_z2){
265 operators_ = &(colSettings->operators);
266 }
267}
268
269template <class Master_matrix>
270template <class Container_type>
271inline Vector_column<Master_matrix>::Vector_column(const Container_type& nonZeroRowIndices,
272 Column_settings* colSettings)
273 : ra_opt(),
274 dim_opt(nonZeroRowIndices.size() == 0 ? 0 : nonZeroRowIndices.size() - 1),
275 chain_opt(),
276 column_(nonZeroRowIndices.size(), nullptr),
277 operators_(nullptr),
278 cellPool_(&(colSettings->cellConstructor))
279{
280 static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type,
281 "Constructor not available for chain columns, please specify the dimension of the chain.");
282
283 if constexpr (!Master_matrix::Option_list::is_z2){
284 operators_ = &(colSettings->operators);
285 }
286
287 index i = 0;
288 if constexpr (Master_matrix::Option_list::is_z2) {
289 for (id_index id : nonZeroRowIndices) {
290 _update_cell(id, i++);
291 }
292 } else {
293 for (const auto& p : nonZeroRowIndices) {
294 _update_cell(operators_->get_value(p.second), p.first, i++);
295 }
296 }
297}
298
299template <class Master_matrix>
300template <class Container_type, class Row_container_type>
301inline Vector_column<Master_matrix>::Vector_column(index columnIndex,
302 const Container_type& nonZeroRowIndices,
303 Row_container_type* rowContainer,
304 Column_settings* colSettings)
305 : ra_opt(columnIndex, rowContainer),
306 dim_opt(nonZeroRowIndices.size() == 0 ? 0 : nonZeroRowIndices.size() - 1),
307 chain_opt([&] {
308 if constexpr (Master_matrix::Option_list::is_z2) {
309 return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : *std::prev(nonZeroRowIndices.end());
310 } else {
311 return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : std::prev(nonZeroRowIndices.end())->first;
312 }
313 }()),
314 column_(nonZeroRowIndices.size(), nullptr),
315 operators_(nullptr),
316 cellPool_(&(colSettings->cellConstructor))
317{
318 static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type,
319 "Constructor not available for chain columns, please specify the dimension of the chain.");
320
321 if constexpr (!Master_matrix::Option_list::is_z2){
322 operators_ = &(colSettings->operators);
323 }
324
325 index i = 0;
326 if constexpr (Master_matrix::Option_list::is_z2) {
327 for (id_index id : nonZeroRowIndices) {
328 _update_cell(id, i++);
329 }
330 } else {
331 for (const auto& p : nonZeroRowIndices) {
332 _update_cell(operators_->get_value(p.second), p.first, i++);
333 }
334 }
335}
336
337template <class Master_matrix>
338template <class Container_type>
339inline Vector_column<Master_matrix>::Vector_column(const Container_type& nonZeroRowIndices,
340 dimension_type dimension,
341 Column_settings* colSettings)
342 : ra_opt(),
343 dim_opt(dimension),
344 chain_opt([&] {
345 if constexpr (Master_matrix::Option_list::is_z2) {
346 return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : *std::prev(nonZeroRowIndices.end());
347 } else {
348 return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : std::prev(nonZeroRowIndices.end())->first;
349 }
350 }()),
351 column_(nonZeroRowIndices.size(), nullptr),
352 operators_(nullptr),
353 cellPool_(&(colSettings->cellConstructor))
354{
355 if constexpr (!Master_matrix::Option_list::is_z2){
356 operators_ = &(colSettings->operators);
357 }
358
359 index i = 0;
360 if constexpr (Master_matrix::Option_list::is_z2) {
361 for (id_index id : nonZeroRowIndices) {
362 _update_cell(id, i++);
363 }
364 } else {
365 for (const auto& p : nonZeroRowIndices) {
366 _update_cell(operators_->get_value(p.second), p.first, i++);
367 }
368 }
369}
370
371template <class Master_matrix>
372template <class Container_type, class Row_container_type>
373inline Vector_column<Master_matrix>::Vector_column(
374 index columnIndex,
375 const Container_type& nonZeroRowIndices,
376 dimension_type dimension,
377 Row_container_type* rowContainer,
378 Column_settings* colSettings)
379 : ra_opt(columnIndex, rowContainer),
380 dim_opt(dimension),
381 chain_opt([&] {
382 if constexpr (Master_matrix::Option_list::is_z2) {
383 return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : *std::prev(nonZeroRowIndices.end());
384 } else {
385 return nonZeroRowIndices.begin() == nonZeroRowIndices.end() ? -1 : std::prev(nonZeroRowIndices.end())->first;
386 }
387 }()),
388 column_(nonZeroRowIndices.size(), nullptr),
389 operators_(nullptr),
390 cellPool_(&(colSettings->cellConstructor))
391{
392 if constexpr (!Master_matrix::Option_list::is_z2){
393 operators_ = &(colSettings->operators);
394 }
395
396 index i = 0;
397 if constexpr (Master_matrix::Option_list::is_z2) {
398 for (id_index id : nonZeroRowIndices) {
399 _update_cell(id, i++);
400 }
401 } else {
402 for (const auto& p : nonZeroRowIndices) {
403 _update_cell(operators_->get_value(p.second), p.first, i++);
404 }
405 }
406}
407
408template <class Master_matrix>
409inline Vector_column<Master_matrix>::Vector_column(const Vector_column& column,
410 Column_settings* colSettings)
411 : ra_opt(),
412 dim_opt(static_cast<const dim_opt&>(column)),
413 chain_opt(static_cast<const chain_opt&>(column)),
414 column_(column.column_.size(), nullptr),
415 erasedValues_(column.erasedValues_),
416 operators_(colSettings == nullptr ? column.operators_ : nullptr),
417 cellPool_(colSettings == nullptr ? column.cellPool_ : &(colSettings->cellConstructor))
418{
419 static_assert(!Master_matrix::Option_list::has_row_access,
420 "Simple copy constructor not available when row access option enabled. Please specify the new column "
421 "index and the row container.");
422
423 if constexpr (!Master_matrix::Option_list::is_z2){
424 if (colSettings != nullptr) operators_ = &(colSettings->operators);
425 }
426
427 index i = 0;
428 for (const Cell* cell : column.column_) {
429 if constexpr (Master_matrix::Option_list::is_z2) {
430 _update_cell(cell->get_row_index(), i++);
431 } else {
432 _update_cell(cell->get_element(), cell->get_row_index(), i++);
433 }
434 }
435}
436
437template <class Master_matrix>
438template <class Row_container_type>
439inline Vector_column<Master_matrix>::Vector_column(const Vector_column& column, index columnIndex,
440 Row_container_type* rowContainer,
441 Column_settings* colSettings)
442 : ra_opt(columnIndex, rowContainer),
443 dim_opt(static_cast<const dim_opt&>(column)),
444 chain_opt(static_cast<const chain_opt&>(column)),
445 column_(column.column_.size(), nullptr),
446 erasedValues_(column.erasedValues_),
447 operators_(colSettings == nullptr ? column.operators_ : nullptr),
448 cellPool_(colSettings == nullptr ? column.cellPool_ : &(colSettings->cellConstructor))
449{
450 if constexpr (!Master_matrix::Option_list::is_z2){
451 if (colSettings != nullptr) operators_ = &(colSettings->operators);
452 }
453
454 index i = 0;
455 for (const Cell* cell : column.column_) {
456 if constexpr (Master_matrix::Option_list::is_z2) {
457 _update_cell(cell->get_row_index(), i++);
458 } else {
459 _update_cell(cell->get_element(), cell->get_row_index(), i++);
460 }
461 }
462}
463
464template <class Master_matrix>
465inline Vector_column<Master_matrix>::Vector_column(Vector_column&& column) noexcept
466 : ra_opt(std::move(static_cast<ra_opt&>(column))),
467 dim_opt(std::move(static_cast<dim_opt&>(column))),
468 chain_opt(std::move(static_cast<chain_opt&>(column))),
469 column_(std::move(column.column_)),
470 erasedValues_(std::move(column.erasedValues_)),
471 operators_(std::exchange(column.operators_, nullptr)),
472 cellPool_(std::exchange(column.cellPool_, nullptr))
473{}
474
475template <class Master_matrix>
476inline Vector_column<Master_matrix>::~Vector_column()
477{
478 for (auto* cell : column_) {
479 _delete_cell(cell);
480 }
481}
482
483template <class Master_matrix>
484inline std::vector<typename Vector_column<Master_matrix>::Field_element_type>
485Vector_column<Master_matrix>::get_content(int columnLength) const
486{
487 if (columnLength < 0 && column_.size() > 0)
488 columnLength = column_.back()->get_row_index() + 1;
489 else if (columnLength < 0)
490 return std::vector<Field_element_type>();
491
492 std::vector<Field_element_type> container(columnLength, 0);
493 for (auto it = column_.begin(); it != column_.end() && (*it)->get_row_index() < static_cast<id_index>(columnLength);
494 ++it) {
495 if constexpr (!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type) {
496 if (erasedValues_.find((*it)->get_row_index()) != erasedValues_.end()) continue;
497 }
498 if constexpr (Master_matrix::Option_list::is_z2) {
499 container[(*it)->get_row_index()] = 1;
500 } else {
501 container[(*it)->get_row_index()] = (*it)->get_element();
502 }
503 }
504 return container;
505}
506
507template <class Master_matrix>
508inline bool Vector_column<Master_matrix>::is_non_zero(id_index rowIndex) const
509{
510 if constexpr (!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type)
511 if (erasedValues_.find(rowIndex) != erasedValues_.end()) return false;
512
513 Cell cell(rowIndex);
514 return std::binary_search(column_.begin(), column_.end(), &cell,
515 [](const Cell* a, const Cell* b) { return a->get_row_index() < b->get_row_index(); });
516}
517
518template <class Master_matrix>
519inline bool Vector_column<Master_matrix>::is_empty() const
520{
521 if constexpr (!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type) {
522 return column_.size() == erasedValues_.size(); // assumes that erasedValues is always a subset of column_, which is
523 // wrong if someone cleared an non exitsing value...
524 } else {
525 return column_.empty();
526 }
527}
528
529template <class Master_matrix>
530inline std::size_t Vector_column<Master_matrix>::size() const
531{
532 if constexpr (!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type) {
533 return column_.size() - erasedValues_.size(); // assumes that erasedValues is always a subset of column_, which is
534 // wrong if someone cleared an non exitsing value...
535 } else {
536 return column_.size();
537 }
538}
539
540template <class Master_matrix>
541template <class Map_type>
542inline void Vector_column<Master_matrix>::reorder(const Map_type& valueMap,
543 [[maybe_unused]] index columnIndex)
544{
545 static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type,
546 "Method not available for chain columns.");
547
548 if (erasedValues_.empty()) { // to avoid useless push_backs.
549 for (Cell* cell : column_) {
550 if constexpr (Master_matrix::Option_list::has_row_access) {
551 ra_opt::unlink(cell);
552 if (columnIndex != static_cast<index>(-1)) cell->set_column_index(columnIndex);
553 }
554 cell->set_row_index(valueMap.at(cell->get_row_index()));
555 if constexpr (Master_matrix::Option_list::has_intrusive_rows && Master_matrix::Option_list::has_row_access)
556 ra_opt::insert_cell(cell->get_row_index(), cell);
557 }
558
559 // all cells have to be deleted first, to avoid problem with insertion when row is a set
560 if constexpr (!Master_matrix::Option_list::has_intrusive_rows && Master_matrix::Option_list::has_row_access) {
561 for (Cell* cell : column_) {
562 ra_opt::insert_cell(cell->get_row_index(), cell);
563 }
564 }
565
566 std::sort(column_.begin(), column_.end(), [](const Cell* c1, const Cell* c2) { return *c1 < *c2; });
567 } else {
568 Column_type newColumn;
569 for (Cell* cell : column_) {
570 if (erasedValues_.find(cell->get_row_index()) == erasedValues_.end()) {
571 if constexpr (Master_matrix::Option_list::has_row_access) {
572 ra_opt::unlink(cell);
573 if (columnIndex != static_cast<index>(-1)) cell->set_column_index(columnIndex);
574 }
575 cell->set_row_index(valueMap.at(cell->get_row_index()));
576 newColumn.push_back(cell);
577 if constexpr (Master_matrix::Option_list::has_intrusive_rows && Master_matrix::Option_list::has_row_access)
578 ra_opt::insert_cell(cell->get_row_index(), cell);
579 } else {
580 _delete_cell(cell);
581 }
582 }
583 // all cells have to be deleted first, to avoid problem with insertion when row is a set
584 if constexpr (!Master_matrix::Option_list::has_intrusive_rows && Master_matrix::Option_list::has_row_access) {
585 for (Cell* cell : column_) {
586 ra_opt::insert_cell(cell->get_row_index(), cell);
587 }
588 }
589 std::sort(newColumn.begin(), newColumn.end(), [](const Cell* c1, const Cell* c2) { return *c1 < *c2; });
590 erasedValues_.clear();
591 column_.swap(newColumn);
592 }
593}
594
595template <class Master_matrix>
596inline void Vector_column<Master_matrix>::clear()
597{
598 static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type,
599 "Method not available for chain columns as a base element should not be empty.");
600
601 for (auto* cell : column_) {
602 if constexpr (Master_matrix::Option_list::has_row_access) ra_opt::unlink(cell);
603 cellPool_->destroy(cell);
604 }
605
606 column_.clear();
607 erasedValues_.clear();
608}
609
610template <class Master_matrix>
611inline void Vector_column<Master_matrix>::clear(id_index rowIndex)
612{
613 static_assert(!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type,
614 "Method not available for chain columns.");
615
616 erasedValues_.insert(rowIndex);
617}
618
619template <class Master_matrix>
620inline typename Vector_column<Master_matrix>::id_index
621Vector_column<Master_matrix>::get_pivot()
622{
623 static_assert(Master_matrix::isNonBasic,
624 "Method not available for base columns."); // could technically be, but is the notion usefull then?
625
626 if constexpr (Master_matrix::Option_list::is_of_boundary_type) {
627 if (column_.empty()) return -1;
628 if (erasedValues_.empty()) return column_.back()->get_row_index();
629
630 auto it = erasedValues_.find(column_.back()->get_row_index());
631 while (!column_.empty() && it != erasedValues_.end()) {
632 erasedValues_.erase(it);
633 _delete_cell(column_.back());
634 column_.pop_back();
635 if (!column_.empty()) it = erasedValues_.find(column_.back()->get_row_index());
636 }
637
638 if (column_.empty()) return -1;
639 return column_.back()->get_row_index();
640 } else {
641 return chain_opt::get_pivot();
642 }
643}
644
645template <class Master_matrix>
646inline typename Vector_column<Master_matrix>::Field_element_type
647Vector_column<Master_matrix>::get_pivot_value()
648{
649 static_assert(Master_matrix::isNonBasic,
650 "Method not available for base columns."); // could technically be, but is the notion usefull then?
651
652 if constexpr (Master_matrix::Option_list::is_z2) {
653 return 1;
654 } else {
655 if constexpr (Master_matrix::Option_list::is_of_boundary_type) {
656 if (column_.empty()) return 0;
657 if (erasedValues_.empty()) return column_.back()->get_element();
658
659 auto it = erasedValues_.find(column_.back()->get_row_index());
660 while (!column_.empty() && it != erasedValues_.end()) {
661 erasedValues_.erase(it);
662 _delete_cell(column_.back());
663 column_.pop_back();
664 if (!column_.empty()) it = erasedValues_.find(column_.back()->get_row_index());
665 }
666
667 if (column_.empty()) return 0;
668 return column_.back()->get_element();
669 } else {
670 if (chain_opt::get_pivot() == static_cast<id_index>(-1)) return Field_element_type();
671 for (const Cell* cell : column_) {
672 if (cell->get_row_index() == chain_opt::get_pivot()) return cell->get_element();
673 }
674 return Field_element_type(); // should never happen if chain column is used properly
675 }
676 }
677}
678
679template <class Master_matrix>
680inline typename Vector_column<Master_matrix>::iterator
681Vector_column<Master_matrix>::begin() noexcept
682{
683 return column_.begin();
684}
685
686template <class Master_matrix>
687inline typename Vector_column<Master_matrix>::const_iterator
688Vector_column<Master_matrix>::begin() const noexcept
689{
690 return column_.begin();
691}
692
693template <class Master_matrix>
694inline typename Vector_column<Master_matrix>::iterator
695Vector_column<Master_matrix>::end() noexcept
696{
697 return column_.end();
698}
699
700template <class Master_matrix>
701inline typename Vector_column<Master_matrix>::const_iterator
702Vector_column<Master_matrix>::end() const noexcept
703{
704 return column_.end();
705}
706
707template <class Master_matrix>
708inline typename Vector_column<Master_matrix>::reverse_iterator
709Vector_column<Master_matrix>::rbegin() noexcept
710{
711 return column_.rbegin();
712}
713
714template <class Master_matrix>
715inline typename Vector_column<Master_matrix>::const_reverse_iterator
716Vector_column<Master_matrix>::rbegin() const noexcept
717{
718 return column_.rbegin();
719}
720
721template <class Master_matrix>
722inline typename Vector_column<Master_matrix>::reverse_iterator
723Vector_column<Master_matrix>::rend() noexcept
724{
725 return column_.rend();
726}
727
728template <class Master_matrix>
729inline typename Vector_column<Master_matrix>::const_reverse_iterator
730Vector_column<Master_matrix>::rend() const noexcept
731{
732 return column_.rend();
733}
734
735template <class Master_matrix>
736template <class Cell_range>
737inline Vector_column<Master_matrix>& Vector_column<Master_matrix>::operator+=(
738 const Cell_range& column)
739{
740 static_assert((!Master_matrix::isNonBasic || std::is_same_v<Cell_range, Vector_column>),
741 "For boundary columns, the range has to be a column of same type to help ensure the validity of the "
742 "base element."); // could be removed, if we give the responsability to the user.
743 static_assert((!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type),
744 "For chain columns, the given column cannot be constant.");
745
746 _add(column);
747
748 return *this;
749}
750
751template <class Master_matrix>
752inline Vector_column<Master_matrix>& Vector_column<Master_matrix>::operator+=(
753 Vector_column& column)
754{
755 if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) {
756 // assumes that the addition never zeros out this column.
757 if (_add(column)) {
758 chain_opt::swap_pivots(column);
759 dim_opt::swap_dimension(column);
760 }
761 } else {
762 _add(column);
763 }
764
765 return *this;
766}
767
768template <class Master_matrix>
769inline Vector_column<Master_matrix>& Vector_column<Master_matrix>::operator*=(
770 unsigned int v)
771{
772 if constexpr (Master_matrix::Option_list::is_z2) {
773 if (v % 2 == 0) {
774 if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) {
775 throw std::invalid_argument("A chain column should not be multiplied by 0.");
776 } else {
777 clear();
778 }
779 }
780 } else {
781 Field_element_type val = operators_->get_value(v);
782
783 if (val == Field_operators::get_additive_identity()) {
784 if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) {
785 throw std::invalid_argument("A chain column should not be multiplied by 0.");
786 } else {
787 clear();
788 }
789 return *this;
790 }
791
792 if (val == Field_operators::get_multiplicative_identity()) return *this;
793
794 for (Cell* cell : column_) {
795 operators_->multiply_inplace(cell->get_element(), val);
796 if constexpr (Master_matrix::Option_list::has_row_access) ra_opt::update_cell(*cell);
797 }
798 }
799
800 return *this;
801}
802
803template <class Master_matrix>
804template <class Cell_range>
805inline Vector_column<Master_matrix>& Vector_column<Master_matrix>::multiply_target_and_add(
806 const Field_element_type& val, const Cell_range& column)
807{
808 static_assert((!Master_matrix::isNonBasic || std::is_same_v<Cell_range, Vector_column>),
809 "For boundary columns, the range has to be a column of same type to help ensure the validity of the "
810 "base element."); // could be removed, if we give the responsability to the user.
811 static_assert((!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type),
812 "For chain columns, the given column cannot be constant.");
813
814 if constexpr (Master_matrix::Option_list::is_z2) {
815 if (val) {
816 _add(column);
817 } else {
818 clear();
819 _add(column);
820 }
821 } else {
822 _multiply_target_and_add(val, column);
823 }
824
825 return *this;
826}
827
828template <class Master_matrix>
829inline Vector_column<Master_matrix>& Vector_column<Master_matrix>::multiply_target_and_add(
830 const Field_element_type& val, Vector_column& column)
831{
832 if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) {
833 // assumes that the addition never zeros out this column.
834 if constexpr (Master_matrix::Option_list::is_z2) {
835 if (val) {
836 if (_add(column)) {
837 chain_opt::swap_pivots(column);
838 dim_opt::swap_dimension(column);
839 }
840 } else {
841 throw std::invalid_argument("A chain column should not be multiplied by 0.");
842 }
843 } else {
844 if (_multiply_target_and_add(val, column)) {
845 chain_opt::swap_pivots(column);
846 dim_opt::swap_dimension(column);
847 }
848 }
849 } else {
850 if constexpr (Master_matrix::Option_list::is_z2) {
851 if (val) {
852 _add(column);
853 } else {
854 clear();
855 _add(column);
856 }
857 } else {
858 _multiply_target_and_add(val, column);
859 }
860 }
861
862 return *this;
863}
864
865template <class Master_matrix>
866template <class Cell_range>
867inline Vector_column<Master_matrix>& Vector_column<Master_matrix>::multiply_source_and_add(
868 const Cell_range& column, const Field_element_type& val)
869{
870 static_assert((!Master_matrix::isNonBasic || std::is_same_v<Cell_range, Vector_column>),
871 "For boundary columns, the range has to be a column of same type to help ensure the validity of the "
872 "base element."); // could be removed, if we give the responsability to the user.
873 static_assert((!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type),
874 "For chain columns, the given column cannot be constant.");
875
876 if constexpr (Master_matrix::Option_list::is_z2) {
877 if (val) {
878 _add(column);
879 }
880 } else {
881 _multiply_source_and_add(column, val);
882 }
883
884 return *this;
885}
886
887template <class Master_matrix>
888inline Vector_column<Master_matrix>& Vector_column<Master_matrix>::multiply_source_and_add(
889 Vector_column& column, const Field_element_type& val)
890{
891 if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) {
892 // assumes that the addition never zeros out this column.
893 if constexpr (Master_matrix::Option_list::is_z2) {
894 if (val) {
895 if (_add(column)) {
896 chain_opt::swap_pivots(column);
897 dim_opt::swap_dimension(column);
898 }
899 }
900 } else {
901 if (_multiply_source_and_add(column, val)) {
902 chain_opt::swap_pivots(column);
903 dim_opt::swap_dimension(column);
904 }
905 }
906 } else {
907 if constexpr (Master_matrix::Option_list::is_z2) {
908 if (val) {
909 _add(column);
910 }
911 } else {
912 _multiply_source_and_add(column, val);
913 }
914 }
915
916 return *this;
917}
918
919template <class Master_matrix>
920inline Vector_column<Master_matrix>& Vector_column<Master_matrix>::operator=(
921 const Vector_column& other)
922{
923 static_assert(!Master_matrix::Option_list::has_row_access, "= assignement not enabled with row access option.");
924
925 dim_opt::operator=(other);
926 chain_opt::operator=(other);
927
928 auto tmpPool = cellPool_;
929 cellPool_ = other.cellPool_;
930
931 while (column_.size() > other.column_.size()) {
932 if (column_.back() != nullptr) {
933 if constexpr (Master_matrix::Option_list::has_row_access) ra_opt::unlink(column_.back());
934 tmpPool->destroy(column_.back());
935 }
936 column_.pop_back();
937 }
938
939 column_.resize(other.column_.size(), nullptr);
940 index i = 0;
941 for (const Cell* cell : other.column_) {
942 if (column_[i] != nullptr) {
943 if constexpr (Master_matrix::Option_list::has_row_access) ra_opt::unlink(column_[i]);
944 tmpPool->destroy(column_[i]);
945 }
946 if constexpr (Master_matrix::Option_list::is_z2) {
947 _update_cell(cell->get_row_index(), i++);
948 } else {
949 _update_cell(cell->get_element(), cell->get_row_index(), i++);
950 }
951 }
952 erasedValues_ = other.erasedValues_;
953 operators_ = other.operators_;
954
955 return *this;
956}
957
958template <class Master_matrix>
959inline std::size_t Vector_column<Master_matrix>::compute_hash_value()
960{
961 std::size_t seed = 0;
962 for (Cell* cell : column_) {
963 if (erasedValues_.find(cell->get_row_index()) == erasedValues_.end()){
964 seed ^= std::hash<unsigned int>()(cell->get_row_index() * static_cast<unsigned int>(cell->get_element())) +
965 0x9e3779b9 + (seed << 6) + (seed >> 2);
966 }
967 }
968 return seed;
969}
970
971template <class Master_matrix>
972inline void Vector_column<Master_matrix>::_delete_cell(Cell* cell)
973{
974 if constexpr (Master_matrix::Option_list::has_row_access) ra_opt::unlink(cell);
975 cellPool_->destroy(cell);
976}
977
978template <class Master_matrix>
979inline void Vector_column<Master_matrix>::_delete_cell(typename Column_type::iterator& it)
980{
981 _delete_cell(*it);
982 ++it;
983}
984
985template <class Master_matrix>
986inline typename Vector_column<Master_matrix>::Cell* Vector_column<Master_matrix>::_insert_cell(
987 const Field_element_type& value, id_index rowIndex, Column_type& column)
988{
989 if constexpr (Master_matrix::Option_list::has_row_access) {
990 Cell* newCell = cellPool_->construct(ra_opt::columnIndex_, rowIndex);
991 newCell->set_element(value);
992 column.push_back(newCell);
993 ra_opt::insert_cell(rowIndex, newCell);
994 return newCell;
995 } else {
996 Cell* newCell = cellPool_->construct(rowIndex);
997 newCell->set_element(value);
998 column.push_back(newCell);
999 return newCell;
1000 }
1001}
1002
1003template <class Master_matrix>
1004inline void Vector_column<Master_matrix>::_insert_cell(id_index rowIndex, Column_type& column)
1005{
1006 if constexpr (Master_matrix::Option_list::has_row_access) {
1007 Cell* newCell = cellPool_->construct(ra_opt::columnIndex_, rowIndex);
1008 column.push_back(newCell);
1009 ra_opt::insert_cell(rowIndex, newCell);
1010 } else {
1011 Cell* newCell = cellPool_->construct(rowIndex);
1012 column.push_back(newCell);
1013 }
1014}
1015
1016template <class Master_matrix>
1017inline void Vector_column<Master_matrix>::_update_cell(const Field_element_type& value,
1018 id_index rowIndex, index position)
1019{
1020 if constexpr (Master_matrix::Option_list::has_row_access) {
1021 Cell* newCell = cellPool_->construct(ra_opt::columnIndex_, rowIndex);
1022 newCell->set_element(value);
1023 column_[position] = newCell;
1024 ra_opt::insert_cell(rowIndex, newCell);
1025 } else {
1026 column_[position] = cellPool_->construct(rowIndex);
1027 column_[position]->set_element(value);
1028 }
1029}
1030
1031template <class Master_matrix>
1032inline void Vector_column<Master_matrix>::_update_cell(id_index rowIndex, index position)
1033{
1034 if constexpr (Master_matrix::Option_list::has_row_access) {
1035 Cell* newCell = cellPool_->construct(ra_opt::columnIndex_, rowIndex);
1036 column_[position] = newCell;
1037 ra_opt::insert_cell(rowIndex, newCell);
1038 } else {
1039 column_[position] = cellPool_->construct(rowIndex);
1040 }
1041}
1042
1043template <class Master_matrix>
1044template <class Cell_range>
1045inline bool Vector_column<Master_matrix>::_add(const Cell_range& column)
1046{
1047 if (column.begin() == column.end()) return false;
1048 if (column_.empty()) { // chain should never enter here.
1049 column_.resize(column.size());
1050 index i = 0;
1051 for (const Cell& cell : column) {
1052 if constexpr (Master_matrix::Option_list::is_z2) {
1053 _update_cell(cell.get_row_index(), i++);
1054 } else {
1055 _update_cell(cell.get_element(), cell.get_row_index(), i++);
1056 }
1057 }
1058 return true;
1059 }
1060
1061 Column_type newColumn;
1062 newColumn.reserve(column_.size() + column.size()); // safe upper bound
1063
1064 auto pivotIsZeroed = _generic_add(
1065 column,
1066 [&](Cell* cellTarget) { newColumn.push_back(cellTarget); },
1067 [&](typename Cell_range::const_iterator& itSource,
1068 [[maybe_unused]] const typename Column_type::iterator& itTarget) {
1069 if constexpr (Master_matrix::Option_list::is_z2) {
1070 _insert_cell(itSource->get_row_index(), newColumn);
1071 } else {
1072 _insert_cell(itSource->get_element(), itSource->get_row_index(), newColumn);
1073 }
1074 },
1075 [&](Field_element_type& targetElement, typename Cell_range::const_iterator& itSource) {
1076 if constexpr (!Master_matrix::Option_list::is_z2)
1077 operators_->add_inplace(targetElement, itSource->get_element());
1078 },
1079 [&](Cell* cellTarget) { newColumn.push_back(cellTarget); });
1080
1081 column_.swap(newColumn);
1082
1083 return pivotIsZeroed;
1084}
1085
1086template <class Master_matrix>
1087template <class Cell_range>
1088inline bool Vector_column<Master_matrix>::_multiply_target_and_add(const Field_element_type& val,
1089 const Cell_range& column)
1090{
1091 if (val == 0u) {
1092 if constexpr (Master_matrix::isNonBasic && !Master_matrix::Option_list::is_of_boundary_type) {
1093 throw std::invalid_argument("A chain column should not be multiplied by 0.");
1094 // this would not only mess up the base, but also the pivots stored.
1095 } else {
1096 clear();
1097 }
1098 }
1099 if (column_.empty()) { // chain should never enter here.
1100 column_.resize(column.size());
1101 index i = 0;
1102 for (const Cell& cell : column) {
1103 if constexpr (Master_matrix::Option_list::is_z2) {
1104 _update_cell(cell.get_row_index(), i++);
1105 } else {
1106 _update_cell(cell.get_element(), cell.get_row_index(), i++);
1107 }
1108 }
1109 if constexpr (std::is_same_v<Cell_range, Vector_column<Master_matrix> >) erasedValues_ = column.erasedValues_;
1110 return true;
1111 }
1112
1113 Column_type newColumn;
1114 newColumn.reserve(column_.size() + column.size()); // safe upper bound
1115
1116 auto pivotIsZeroed = _generic_add(
1117 column,
1118 [&](Cell* cellTarget) {
1119 operators_->multiply_inplace(cellTarget->get_element(), val);
1120 if constexpr (Master_matrix::Option_list::has_row_access) ra_opt::update_cell(*cellTarget);
1121 newColumn.push_back(cellTarget);
1122 },
1123 [&](typename Cell_range::const_iterator& itSource, const typename Column_type::iterator& itTarget) {
1124 _insert_cell(itSource->get_element(), itSource->get_row_index(), newColumn);
1125 },
1126 [&](Field_element_type& targetElement, typename Cell_range::const_iterator& itSource) {
1127 operators_->multiply_and_add_inplace_front(targetElement, val, itSource->get_element());
1128 },
1129 [&](Cell* cellTarget) { newColumn.push_back(cellTarget); });
1130
1131 column_.swap(newColumn);
1132
1133 return pivotIsZeroed;
1134}
1135
1136template <class Master_matrix>
1137template <class Cell_range>
1138inline bool Vector_column<Master_matrix>::_multiply_source_and_add(const Cell_range& column,
1139 const Field_element_type& val)
1140{
1141 if (val == 0u || column.begin() == column.end()) {
1142 return false;
1143 }
1144
1145 Column_type newColumn;
1146 newColumn.reserve(column_.size() + column.size()); // safe upper bound
1147
1148 auto pivotIsZeroed = _generic_add(
1149 column,
1150 [&](Cell* cellTarget) { newColumn.push_back(cellTarget); },
1151 [&](typename Cell_range::const_iterator& itSource, const typename Column_type::iterator& itTarget) {
1152 Cell* newCell = _insert_cell(itSource->get_element(), itSource->get_row_index(), newColumn);
1153 operators_->multiply_inplace(newCell->get_element(), val);
1154 if constexpr (Master_matrix::Option_list::has_row_access) ra_opt::update_cell(*newCell);
1155 },
1156 [&](Field_element_type& targetElement, typename Cell_range::const_iterator& itSource) {
1157 operators_->multiply_and_add_inplace_back(itSource->get_element(), val, targetElement);
1158 },
1159 [&](Cell* cellTarget) { newColumn.push_back(cellTarget); });
1160
1161 column_.swap(newColumn);
1162
1163 return pivotIsZeroed;
1164}
1165
1166template <class Master_matrix>
1167template <class Cell_range, typename F1, typename F2, typename F3, typename F4>
1168inline bool Vector_column<Master_matrix>::_generic_add(const Cell_range& column,
1169 F1&& process_target,
1170 F2&& process_source,
1171 F3&& update_target1,
1172 F4&& update_target2)
1173{
1174 auto updateTargetIterator = [&](typename Column_type::iterator& itTarget){
1175 if constexpr (!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type){
1176 while (itTarget != column_.end() && erasedValues_.find((*itTarget)->get_row_index()) != erasedValues_.end()) {
1177 _delete_cell(*itTarget);
1178 ++itTarget;
1179 }
1180 }
1181 };
1182 auto updateSourceIterator = [&](typename Cell_range::const_iterator& itSource){
1183 if constexpr (std::is_same_v<Cell_range, Vector_column<Master_matrix> > &&
1184 (!Master_matrix::isNonBasic || Master_matrix::Option_list::is_of_boundary_type)) {
1185 while (itSource != column.end() &&
1186 column.erasedValues_.find(itSource->get_row_index()) != column.erasedValues_.end())
1187 ++itSource;
1188 }
1189 };
1190
1191 bool pivotIsZeroed = false;
1192
1193 auto itTarget = column_.begin();
1194 auto itSource = column.begin();
1195 while (itTarget != column_.end() && itSource != column.end()) {
1196 updateTargetIterator(itTarget);
1197 updateSourceIterator(itSource);
1198 if (itTarget == column_.end() || itSource == column.end()) break;
1199
1200 _generic_merge_cell_to_column(*this,
1201 itSource, itTarget,
1202 process_target, process_source, update_target1, update_target2,
1203 pivotIsZeroed);
1204 }
1205
1206 while (itSource != column.end()) {
1207 updateSourceIterator(itSource);
1208 if (itSource == column.end()) break;
1209
1210 process_source(itSource, column_.end());
1211 ++itSource;
1212 }
1213
1214 while (itTarget != column_.end()) {
1215 updateTargetIterator(itTarget);
1216 if (itTarget == column_.end()) break;
1217
1218 process_target(*itTarget);
1219 ++itTarget;
1220 }
1221
1222 erasedValues_.clear();
1223
1224 return pivotIsZeroed;
1225}
1226
1227} // namespace persistence_matrix
1228} // namespace Gudhi
1229
1238template <class Master_matrix>
1239struct std::hash<Gudhi::persistence_matrix::Vector_column<Master_matrix> >
1240{
1241 std::size_t operator()(const Gudhi::persistence_matrix::Vector_column<Master_matrix>& column) const {
1242 return column.compute_hash_value();
1243 }
1244};
1245
1246#endif // PM_VECTOR_COLUMN_H
Column class following the PersistenceMatrixColumn concept.
Definition vector_column.h:53
Contains helper methods for column addition and column hasher.
Gudhi namespace.
Definition SimplicialComplexForAlpha.h:14