73 typedef Persistent_cohomology_column<Simplex_key, Arith_element> Column;
75 typedef typename Column::Cell Cell;
77 typedef boost::intrusive::list<Cell,
78 boost::intrusive::constant_time_size<false>,
79 boost::intrusive::base_hook<base_hook_cam_h> > Hcell;
81 typedef boost::intrusive::set<Column,
82 boost::intrusive::constant_time_size<false> > Cam;
84 typedef std::vector<std::pair<Simplex_key, Arith_element> > A_ds_type;
99 dim_max_(cpx.dimension()),
101 num_simplices_(cpx_->num_simplices()),
102 ds_rank_(num_simplices_),
103 ds_parent_(num_simplices_),
104 ds_repr_(num_simplices_, NULL),
105 dsets_(ds_rank_.data(), ds_parent_.data()),
110 interval_length_policy(&cpx, 0),
113 if (num_simplices_ > std::numeric_limits<Simplex_key>::max()) {
115 throw std::out_of_range(
"The number of simplices is more than Simplex_key type numeric limit.");
117 if (persistence_dim_max) {
124 for (
auto & transverse_ref : transverse_idx_) {
126 transverse_ref.second.row_->clear_and_dispose([&](Cell*p){p->~Cell();});
127 delete transverse_ref.second.row_;
132 struct length_interval {
135 min_length_(min_length) {
138 bool operator()(Simplex_handle sh1, Simplex_handle sh2) {
139 return cpx_->filtration(sh2) - cpx_->filtration(sh1) > min_length_;
143 min_length_ = new_length;
153 coeff_field_.init(charac);
157 coeff_field_.init(charac_min, charac_max);
169 interval_length_policy.set_length(min_interval_length);
171 std::vector<Simplex_key> vertices;
175 dsets_.make_set(cpx_->
key(sh));
177 switch (dim_simplex) {
179 vertices.push_back(idx_fil);
182 update_cohomology_groups_edge(sh);
185 update_cohomology_groups(sh, dim_simplex);
191 if (ds_parent_[key] == key
192 && zero_cocycles_.find(key) == zero_cocycles_.end()) {
193 persistent_pairs_.emplace_back(
197 for (
auto zero_idx : zero_cocycles_) {
198 persistent_pairs_.emplace_back(
202 for (
auto cocycle : transverse_idx_) {
203 persistent_pairs_.emplace_back(
213 void update_cohomology_groups_edge(Simplex_handle sigma) {
215 boost::tie(u, v) = cpx_->endpoints(sigma);
225 auto map_it_u = zero_cocycles_.find(ku);
227 if (map_it_u == zero_cocycles_.end()) {
230 idx_coc_u = map_it_u->second;
233 auto map_it_v = zero_cocycles_.find(kv);
235 if (map_it_v == zero_cocycles_.end()) {
238 idx_coc_v = map_it_v->second;
243 if (interval_length_policy(cpx_->
simplex(idx_coc_v), sigma)) {
244 persistent_pairs_.emplace_back(
248 if (kv != idx_coc_v) {
249 zero_cocycles_.erase(map_it_v);
251 if (kv == dsets_.find_set(kv)) {
252 if (ku != idx_coc_u) {
253 zero_cocycles_.erase(map_it_u);
255 zero_cocycles_[kv] = idx_coc_u;
258 if (interval_length_policy(cpx_->
simplex(idx_coc_u), sigma)) {
259 persistent_pairs_.emplace_back(
263 if (ku != idx_coc_u) {
264 zero_cocycles_.erase(map_it_u);
266 if (ku == dsets_.find_set(ku)) {
267 if (kv != idx_coc_v) {
268 zero_cocycles_.erase(map_it_v);
270 zero_cocycles_[ku] = idx_coc_v;
274 }
else if (dim_max_ > 1) {
282 void annotation_of_the_boundary(
283 std::map<Simplex_key, Arith_element> & map_a_ds, Simplex_handle sigma,
288 typedef std::pair<Column *, int> annotation_t;
289 thread_local std::vector<annotation_t> annotations_in_boundary;
290 annotations_in_boundary.clear();
291 int sign = 1 - 2 * (dim_sigma % 2);
296 for (
auto sh : cpx_->boundary_simplex_range(sigma)) {
300 curr_col = ds_repr_[dsets_.find_set(key)];
301 if (curr_col != NULL) {
302 annotations_in_boundary.emplace_back(curr_col, sign);
308 std::sort(annotations_in_boundary.begin(), annotations_in_boundary.end(),
309 [](annotation_t
const& a, annotation_t
const& b) { return a.first < b.first; });
313 std::pair<typename std::map<Simplex_key, Arith_element>::iterator,
bool> result_insert_a_ds;
315 for (
auto ann_it = annotations_in_boundary.begin(); ann_it != annotations_in_boundary.end(); ) {
316 Column* col = ann_it->first;
317 int mult = ann_it->second;
318 while (++ann_it != annotations_in_boundary.end() && ann_it->first == col) {
319 mult += ann_it->second;
323 for (
auto cell_ref : col->col_) {
324 Arith_element w_y = coeff_field_.times(cell_ref.coefficient_, mult);
327 result_insert_a_ds = map_a_ds.insert(std::pair<Simplex_key, Arith_element>(cell_ref.key_, w_y));
328 if (!(result_insert_a_ds.second)) {
329 result_insert_a_ds.first->second = coeff_field_.
plus_equal(result_insert_a_ds.first->second, w_y);
331 map_a_ds.erase(result_insert_a_ds.first);
343 void update_cohomology_groups(Simplex_handle sigma,
int dim_sigma) {
345 std::map<Simplex_key, Arith_element> map_a_ds;
346 annotation_of_the_boundary(map_a_ds, sigma, dim_sigma);
348 if (map_a_ds.empty()) {
349 if (dim_sigma < dim_max_) {
356 for (
auto map_a_ds_ref : map_a_ds) {
358 std::pair<Simplex_key, Arith_element>(map_a_ds_ref.first,
359 map_a_ds_ref.second));
364 for (
auto a_ds_rit = a_ds.rbegin();
365 (a_ds_rit != a_ds.rend())
367 std::tie(inv_x, charac) = coeff_field_.inverse(a_ds_rit->second, prod);
370 destroy_cocycle(sigma, a_ds, a_ds_rit->first, inv_x, charac);
375 && dim_sigma < dim_max_) {
392 Column * new_col = column_pool_.construct(key);
393 Cell * new_cell = cell_pool_.construct(key, x, new_col);
394 new_col->col_.push_back(*new_cell);
398 cam_.insert(cam_.end(), *new_col);
400 Hcell * new_hcell =
new Hcell;
401 new_hcell->push_back(*new_cell);
402 transverse_idx_[key] = cocycle(charac, new_hcell);
403 ds_repr_[key] = new_col;
414 void destroy_cocycle(Simplex_handle sigma, A_ds_type
const& a_ds,
418 if (interval_length_policy(cpx_->
simplex(death_key), sigma)) {
419 persistent_pairs_.emplace_back(cpx_->
simplex(death_key)
424 auto death_key_row = transverse_idx_.find(death_key);
425 std::pair<typename Cam::iterator, bool> result_insert_cam;
427 auto row_cell_it = death_key_row->second.row_->begin();
429 while (row_cell_it != death_key_row->second.row_->end()) {
431 Arith_element w = coeff_field_.times_minus(inv_x, row_cell_it->coefficient_);
434 Column * curr_col = row_cell_it->self_col_;
437 for (
auto& col_cell : curr_col->col_) {
438 col_cell.base_hook_cam_h::unlink();
442 cam_.erase(cam_.iterator_to(*curr_col));
444 plus_equal_column(*curr_col, a_ds, w);
446 if (curr_col->col_.empty()) {
447 ds_repr_[curr_col->class_key_] = NULL;
448 column_pool_.destroy(curr_col);
451 result_insert_cam = cam_.insert(*curr_col);
452 if (result_insert_cam.second) {
453 for (
auto& col_cell : curr_col->col_) {
455 transverse_idx_[col_cell.key_].row_->push_front(col_cell);
459 dsets_.link(curr_col->class_key_,
460 result_insert_cam.first->class_key_);
462 Simplex_key key_tmp = dsets_.find_set(curr_col->class_key_);
463 ds_repr_[key_tmp] = &(*(result_insert_cam.first));
464 result_insert_cam.first->class_key_ = key_tmp;
466 curr_col->col_.clear_and_dispose([&](Cell*p){cell_pool_.destroy(p);});
467 column_pool_.destroy(curr_col);
479 if (death_key_row->second.characteristics_ == charac) {
480 delete death_key_row->second.row_;
481 transverse_idx_.erase(death_key_row);
483 death_key_row->second.characteristics_ /= charac;
490 void plus_equal_column(Column & target, A_ds_type
const& other
492 auto target_it = target.col_.begin();
493 auto other_it = other.begin();
494 while (target_it != target.col_.end() && other_it != other.end()) {
495 if (target_it->key_ < other_it->first) {
498 if (target_it->key_ > other_it->first) {
499 Cell * cell_tmp = cell_pool_.construct(Cell(other_it->first
502 cell_tmp->coefficient_ = coeff_field_.plus_times_equal(cell_tmp->coefficient_, other_it->second, w);
504 target.col_.insert(target_it, *cell_tmp);
509 target_it->coefficient_ = coeff_field_.plus_times_equal(target_it->coefficient_, other_it->second, w);
511 auto tmp_it = target_it;
514 Cell * tmp_cell_ptr = &(*tmp_it);
515 target.col_.erase(tmp_it);
517 cell_pool_.destroy(tmp_cell_ptr);
525 while (other_it != other.end()) {
526 Cell * cell_tmp = cell_pool_.construct(Cell(other_it->first, coeff_field_.
additive_identity(), &target));
527 cell_tmp->coefficient_ = coeff_field_.plus_times_equal(cell_tmp->coefficient_, other_it->second, w);
528 target.col_.insert(target.col_.end(), *cell_tmp);
537 struct cmp_intervals_by_length {
542 return (sc_->filtration(get < 1 > (p1)) - sc_->filtration(get < 0 > (p1))
543 > sc_->filtration(get < 1 > (p2)) - sc_->filtration(get < 0 > (p2)));
560 cmp_intervals_by_length cmp(cpx_);
561 std::sort(std::begin(persistent_pairs_), std::end(persistent_pairs_), cmp);
562 for (
auto pair : persistent_pairs_) {
563 ostream << get<2>(pair) <<
" " << cpx_->
dimension(get<0>(pair)) <<
" "
565 << cpx_->
filtration(get<1>(pair)) <<
" " << std::endl;
569 void write_output_diagram(std::string diagram_name) {
570 std::ofstream diagram_out(diagram_name.c_str());
571 diagram_out.exceptions(diagram_out.failbit);
572 cmp_intervals_by_length cmp(cpx_);
573 std::sort(std::begin(persistent_pairs_), std::end(persistent_pairs_), cmp);
574 for (
auto pair : persistent_pairs_) {
575 diagram_out << cpx_->
dimension(get<0>(pair)) <<
" "
577 << cpx_->
filtration(get<1>(pair)) << std::endl;
589 for (
auto pair : persistent_pairs_) {
607 for (
auto pair : persistent_pairs_) {
610 if (cpx_->
dimension(get<0>(pair)) == dimension) {
628 for (
auto pair : persistent_pairs_) {
650 for (
auto pair : persistent_pairs_) {
656 if (cpx_->
dimension(get<0>(pair)) == dimension) {
669 return persistent_pairs_;
676 std::vector< std::pair< Filtration_value , Filtration_value > >
678 std::vector< std::pair< Filtration_value , Filtration_value > > result;
680 for (
auto && pair : persistent_pairs_) {
681 if (cpx_->
dimension(get<0>(pair)) == dimension) {
699 characteristics_(characteristics) {
710 size_t num_simplices_;
716 std::vector<int> ds_rank_;
717 std::vector<Simplex_key> ds_parent_;
718 std::vector<Column *> ds_repr_;
719 boost::disjoint_sets<int *, Simplex_key *> dsets_;
725 std::unordered_map<Simplex_key, Simplex_key> zero_cocycles_;
727 std::map<Simplex_key, cocycle> transverse_idx_;
729 std::vector<Persistent_interval> persistent_pairs_;
730 length_interval interval_length_policy;
732 Simple_object_pool<Column> column_pool_;
733 Simple_object_pool<Cell> cell_pool_;