diff options
Diffstat (limited to 'matching/include/bifiltration.hpp')
| -rw-r--r-- | matching/include/bifiltration.hpp | 421 |
1 files changed, 421 insertions, 0 deletions
diff --git a/matching/include/bifiltration.hpp b/matching/include/bifiltration.hpp new file mode 100644 index 0000000..9e2a82e --- /dev/null +++ b/matching/include/bifiltration.hpp @@ -0,0 +1,421 @@ +namespace md { + + template<class Real> + void Bifiltration<Real>::init() + { + auto lower_left = max_point<Real>(); + auto upper_right = min_point<Real>(); + 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<Real>(lower_left, upper_right); + } + + template<class Real> + Bifiltration<Real>::Bifiltration(const std::string& fname) + { + 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); + } + + BifiltrationFormat input_format; + + std::string s; + + while(ignore_line(s)) { + std::getline(ifstr, s); + } + + if (s == "bifiltration") { + input_format = BifiltrationFormat::rivet; + } else if (s == "bifiltration_phat_like") { + input_format = BifiltrationFormat::phat_like; + } else { + std::cerr << "Unknown format: '" << s << "' in file " << fname << std::endl; + throw std::runtime_error("unknown bifiltration format"); + } + + switch(input_format) { + case BifiltrationFormat::rivet : + rivet_format_reader(ifstr); + break; + case BifiltrationFormat::phat_like : + phat_like_format_reader(ifstr); + break; + } + + ifstr.close(); + + init(); + } + + template<class Real> + void Bifiltration<Real>::rivet_format_reader(std::ifstream& ifstr) + { + std::string s; + // read axes names, ignore them + std::getline(ifstr, s); + std::getline(ifstr, s); + + Index index = 0; + while(std::getline(ifstr, s)) { + if (!ignore_line(s)) { + simplices_.emplace_back(index++, s, BifiltrationFormat::rivet); + } + } + } + + template<class Real> + void Bifiltration<Real>::phat_like_format_reader(std::ifstream& ifstr) + { + spd::debug("Enter phat_like_format_reader"); + // read stream line by line; do not use >> operator + std::string s; + std::getline(ifstr, s); + + // first line contains number of simplices + long n_simplices = std::stol(s); + + // all other lines represent a simplex + Index index = 0; + while(index < n_simplices) { + std::getline(ifstr, s); + if (!ignore_line(s)) { + simplices_.emplace_back(index++, s, BifiltrationFormat::phat_like); + } + } + spd::debug("Read {} simplices from file", n_simplices); + } + + template<class Real> + void Bifiltration<Real>::scale(Real lambda) + { + for(auto& s : simplices_) { + s.scale(lambda); + } + init(); + } + + template<class Real> + void Bifiltration<Real>::sanity_check() const + { +#ifdef DEBUG + spd::debug("Enter Bifiltration<Real>::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<Real>::sanity_check"); +#endif + } + + template<class Real> + Diagram<Real> Bifiltration<Real>::weighted_slice_diagram(const DualPoint<Real>& line, int dim) const + { + DiagramKeeper<Real> dgm; + + // make a copy for now; I want slice_diagram to be const + std::vector<Simplex<Real>> 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<Real>& a, const Simplex<Real>& 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(); + constexpr Real real_inf = std::numeric_limits<Real>::infinity(); + 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); + bool is_finite_pair = new_pair.second != phat::k_infinity_index; + Real birth = simplices.at(new_pair.first).value(); + Real death = is_finite_pair ? simplices.at(new_pair.second).value() : real_inf; + 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.get_diagram(dim); + } + + template<class Real> + Box<Real> Bifiltration<Real>::bounding_box() const + { + return bounding_box_; + } + + template<class Real> + Real Bifiltration<Real>::minimal_coordinate() const + { + return std::min(bounding_box_.lower_left().x, bounding_box_.lower_left().y); + } + + template<class Real> + void Bifiltration<Real>::translate(Real a) + { + bounding_box_.translate(a); + for(auto& simplex : simplices_) { + simplex.translate(a); + } + } + + template<class Real> + Real Bifiltration<Real>::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; + } + + template<class Real> + Real Bifiltration<Real>::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; + } + + template<class Real> + Real Bifiltration<Real>::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; + } + + template<class Real> + Real Bifiltration<Real>::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; + } + + template<class Real> + void Bifiltration<Real>::add_simplex(Index _id, Point<Real> birth, int _dim, const Column& _bdry) + { + simplices_.emplace_back(_id, birth, _dim, _bdry); + } + + template<class Real> + void Bifiltration<Real>::save(const std::string& filename, md::BifiltrationFormat format) + { + switch(format) { + case BifiltrationFormat::rivet: + throw std::runtime_error("Not implemented"); + break; + case BifiltrationFormat::phat_like: { + 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; + } + } + + template<class Real> + void Bifiltration<Real>::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 + } + + template<class Real> + std::ostream& operator<<(std::ostream& os, const Bifiltration<Real>& bif) + { + os << "Bifiltration [" << std::endl; + for(const auto& s : bif.simplices()) { + os << s << std::endl; + } + os << "]" << std::endl; + return os; + } + + template<class Real> + BifiltrationProxy<Real>::BifiltrationProxy(const Bifiltration<Real>& bif, int dim) + : + dim_(dim), + bif_(bif) + { + cache_positions(); + } + + template<class Real> + void BifiltrationProxy<Real>::cache_positions() const + { + cached_positions_.clear(); + for(const auto& simplex : bif_.simplices()) { + if (simplex.dim() == dim_ or simplex.dim() == dim_ + 1) + cached_positions_.push_back(simplex.position()); + } + } + + template<class Real> + PointVec<Real> + BifiltrationProxy<Real>::positions() const + { + if (cached_positions_.empty()) { + cache_positions(); + } + return cached_positions_; + } + + // translate all points by vector (a,a) + template<class Real> + void BifiltrationProxy<Real>::translate(Real a) + { + bif_.translate(a); + } + + // return minimal value of x- and y-coordinates + // among all simplices + template<class Real> + Real BifiltrationProxy<Real>::minimal_coordinate() const + { + return bif_.minimal_coordinate(); + } + + // return box that contains positions of all simplices + template<class Real> + Box<Real> BifiltrationProxy<Real>::bounding_box() const + { + return bif_.bounding_box(); + } + + template<class Real> + Real BifiltrationProxy<Real>::max_x() const + { + return bif_.max_x(); + } + + template<class Real> + Real BifiltrationProxy<Real>::max_y() const + { + return bif_.max_y(); + } + + template<class Real> + Real BifiltrationProxy<Real>::min_x() const + { + return bif_.min_x(); + } + + template<class Real> + Real BifiltrationProxy<Real>::min_y() const + { + return bif_.min_y(); + } + + + template<class Real> + Diagram<Real> BifiltrationProxy<Real>::weighted_slice_diagram(const DualPoint<Real>& slice) const + { + return bif_.weighted_slice_diagram(slice, dim_); + } + +} + |