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