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