diff options
Diffstat (limited to 'matching/src/bifiltration.cpp')
| -rw-r--r-- | matching/src/bifiltration.cpp | 309 |
1 files changed, 309 insertions, 0 deletions
diff --git a/matching/src/bifiltration.cpp b/matching/src/bifiltration.cpp new file mode 100644 index 0000000..429a9a8 --- /dev/null +++ b/matching/src/bifiltration.cpp @@ -0,0 +1,309 @@ +#include <iostream> +#include <fstream> +#include <sstream> +#include <cassert> + +#include<phat/boundary_matrix.h> +#include<phat/compute_persistence_pairs.h> + +#include "spdlog/spdlog.h" +#include "spdlog/fmt/fmt.h" +#include "spdlog/fmt/ostr.h" + +#include "common_util.h" +#include "bifiltration.h" + +namespace spd = spdlog; + +namespace md { + + void Bifiltration::init() + { + Point lower_left = max_point(); + Point upper_right = min_point(); + for(const auto& simplex : simplices_) { + lower_left = greatest_lower_bound(lower_left, simplex.position()); + upper_right = least_upper_bound(upper_right, simplex.position()); + maximal_dim_ = std::max(maximal_dim_, simplex.dim()); + } + bounding_box_ = Box(lower_left, upper_right); + } + + Bifiltration::Bifiltration(const std::string& fname, BifiltrationFormat input_format) + { + std::ifstream ifstr{fname.c_str()}; + if (!ifstr.good()) { + std::string error_message = fmt::format("Cannot open file {0}", fname); + std::cerr << error_message << std::endl; + throw std::runtime_error(error_message); + } + + switch(input_format) { + case BifiltrationFormat::rivet : + rivet_format_reader(ifstr); + break; + case BifiltrationFormat::rene : + rene_format_reader(ifstr); + break; + } + init(); + } + + void Bifiltration::rivet_format_reader(std::ifstream& ifstr) + { + std::string s; + std::getline(ifstr, s); + assert(s == std::string("bifiltration")); + std::getline(ifstr, parameter_1_name_); + std::getline(ifstr, parameter_2_name_); + + Index index = 0; + while(std::getline(ifstr, s)) { + if (not ignore_line(s)) + simplices_.emplace_back(index++, s, BifiltrationFormat::rivet); + } + } + + void Bifiltration::rene_format_reader(std::ifstream& ifstr) + { + spd::debug("Enter rene_format_reader"); + // read stream line by line; do not use >> operator + std::string s; + std::getline(ifstr, s); + long n_simplices = std::stol(s); + + for(Index index = 0; index < n_simplices; index++) { + std::getline(ifstr, s); + simplices_.emplace_back(index, s, BifiltrationFormat::rene); + } + spd::debug("Read {} simplices from file", n_simplices); + } + + void Bifiltration::scale(Real lambda) + { + for(auto& s : simplices_) { + s.scale(lambda); + } + init(); + } + + void Bifiltration::sanity_check() const + { +#ifdef DEBUG + spd::debug("Enter Bifiltration::sanity_check"); + // check that boundary has correct number of simplices, + // each bounding simplex has correct dim + // and appears in the filtration before the simplex it bounds + for(const auto& s : simplices_) { + assert(s.dim() >= 0); + assert(s.dim() == 0 or s.dim() + 1 == (int) s.boundary().size()); + for(auto bdry_idx : s.boundary()) { + Simplex bdry_simplex = simplices()[bdry_idx]; + assert(bdry_simplex.dim() == s.dim() - 1); + assert(bdry_simplex.position().is_less(s.position(), false)); + } + } + spd::debug("Exit Bifiltration::sanity_check"); +#endif + } + + DiagramKeeper Bifiltration::weighted_slice_diagram(const DualPoint& line, int /*dim*/) const + { + DiagramKeeper dgm; + + // make a copy for now; I want slice_diagram to be const + std::vector<Simplex> simplices(simplices_); + +// std::vector<Simplex> simplices; +// simplices.reserve(simplices_.size() / 2); +// for(const auto& s : simplices_) { +// if (s.dim() <= dim + 1 and s.dim() >= dim) +// simplices.emplace_back(s); +// } + + spd::debug("Enter slice diagram, line = {}, simplices.size = {}", line, simplices.size()); + + for(auto& simplex : simplices) { + Real value = line.weighted_push(simplex.position()); +// spd::debug("in slice_diagram, simplex = {}, value = {}\n", simplex, value); + simplex.set_value(value); + } + + std::sort(simplices.begin(), simplices.end(), + [](const Simplex& a, const Simplex& b) { return a.value() < b.value(); }); + std::map<Index, Index> index_map; + for(Index i = 0; i < (int) simplices.size(); i++) { + index_map[simplices[i].id()] = i; + } + + phat::boundary_matrix<> phat_matrix; + phat_matrix.set_num_cols(simplices.size()); + std::vector<Index> bd_in_slice_filtration; + for(Index i = 0; i < (int) simplices.size(); i++) { + phat_matrix.set_dim(i, simplices[i].dim()); + bd_in_slice_filtration.clear(); + //std::cout << "new col" << i << std::endl; + for(int j = 0; j < (int) simplices[i].boundary().size(); j++) { + // F[i] contains the indices of its facet wrt to the + // original filtration. We have to express it, however, + // wrt to the filtration along the slice. That is why + // we need the index_map + //std::cout << "Found " << F[i].bd[j] << ", returning " << index_map[F[i].bd[j]] << std::endl; + bd_in_slice_filtration.push_back(index_map[simplices[i].boundary()[j]]); + } + std::sort(bd_in_slice_filtration.begin(), bd_in_slice_filtration.end()); + phat_matrix.set_col(i, bd_in_slice_filtration); + } + phat::persistence_pairs phat_persistence_pairs; + phat::compute_persistence_pairs<phat::twist_reduction>(phat_persistence_pairs, phat_matrix); + + dgm.clear(); + for(long i = 0; i < (long) phat_persistence_pairs.get_num_pairs(); i++) { + std::pair<phat::index, phat::index> new_pair = phat_persistence_pairs.get_pair(i); + Real birth = simplices.at(new_pair.first).value(); + Real death = simplices.at(new_pair.second).value(); + int dim = simplices[new_pair.first].dim(); + assert(dim + 1 == simplices[new_pair.second].dim()); + if (birth != death) { + dgm.add_point(dim, birth, death); + } + } + + spdlog::debug("Exiting slice_diagram, #dgm[0] = {}", dgm.get_diagram(0).size()); + + return dgm; + } + + Box Bifiltration::bounding_box() const + { + return bounding_box_; + } + + Real Bifiltration::minimal_coordinate() const + { + return std::min(bounding_box_.lower_left().x, bounding_box_.lower_left().y); + } + + void Bifiltration::translate(Real a) + { + bounding_box_.translate(a); + for(auto& simplex : simplices_) { + simplex.translate(a); + } + } + + Real Bifiltration::max_x() const + { + if (simplices_.empty()) + return 1; + auto me = std::max_element(simplices_.cbegin(), simplices_.cend(), + [](const auto& s_a, const auto& s_b) { return s_a.position().x < s_b.position().x; }); + assert(me != simplices_.cend()); + return me->position().x; + } + + Real Bifiltration::max_y() const + { + if (simplices_.empty()) + return 1; + auto me = std::max_element(simplices_.cbegin(), simplices_.cend(), + [](const auto& s_a, const auto& s_b) { return s_a.position().y < s_b.position().y; }); + assert(me != simplices_.cend()); + return me->position().y; + } + + Real Bifiltration::min_x() const + { + if (simplices_.empty()) + return 0; + auto me = std::min_element(simplices_.cbegin(), simplices_.cend(), + [](const auto& s_a, const auto& s_b) { return s_a.position().x < s_b.position().x; }); + assert(me != simplices_.cend()); + return me->position().x; + } + + Real Bifiltration::min_y() const + { + if (simplices_.empty()) + return 0; + auto me = std::min_element(simplices_.cbegin(), simplices_.cend(), + [](const auto& s_a, const auto& s_b) { return s_a.position().y < s_b.position().y; }); + assert(me != simplices_.cend()); + return me->position().y; + } + + void Bifiltration::add_simplex(md::Index _id, md::Point birth, int _dim, const md::Column& _bdry) + { + simplices_.emplace_back(_id, birth, _dim, _bdry); + } + + void Bifiltration::save(const std::string& filename, md::BifiltrationFormat format) + { + switch(format) { + case BifiltrationFormat::rivet: + throw std::runtime_error("Not implemented"); + break; + case BifiltrationFormat::rene: { + std::ofstream f(filename); + if (not f.good()) { + std::cerr << "Bifiltration::save: cannot open file " << filename << std::endl; + throw std::runtime_error("Cannot open file for writing "); + } + f << simplices_.size() << "\n"; + + for(const auto& s : simplices_) { + f << s.dim() << " " << s.position().x << " " << s.position().y << " "; + for(int b : s.boundary()) { + f << b << " "; + } + f << std::endl; + } + + } + break; + } + } + + + void Bifiltration::postprocess_rivet_format() + { + std::map<Column, Index> facets_to_ids; + + // fill the map + for(Index i = 0; i < (Index)simplices_.size(); ++i) { + assert(simplices_[i].id() == i); + facets_to_ids[simplices_[i].vertices_] = i; + } + +// for(const auto& s : simplices_) { +// facets_to_ids[s] = s.id(); +// } + + // main loop + for(auto& s : simplices_) { + assert(not s.vertices_.empty()); + assert(s.facet_indices_.empty()); + Column facet_indices; + for(Index i = 0; i <= s.dim(); ++i) { + Column facet; + for(Index j : s.vertices_) { + if (j != i) + facet.push_back(j); + } + auto facet_index = facets_to_ids.at(facet); + facet_indices.push_back(facet_index); + } + s.facet_indices_ = facet_indices; + } // loop over simplices + } + + std::ostream& operator<<(std::ostream& os, const Bifiltration& bif) + { + os << "Bifiltration, axes = " << bif.parameter_1_name_ << ", " << bif.parameter_2_name_ << std::endl; + for(const auto& s : bif.simplices()) { + os << s << std::endl; + } + return os; + } +} + |