diff options
Diffstat (limited to 'matching/include')
| -rw-r--r-- | matching/include/bifiltration.h | 31 | ||||
| -rw-r--r-- | matching/include/common_util.h | 43 | ||||
| -rw-r--r-- | matching/include/matching_distance.h | 59 | ||||
| -rw-r--r-- | matching/include/matching_distance.hpp | 784 | ||||
| -rw-r--r-- | matching/include/persistence_module.h | 63 | ||||
| -rw-r--r-- | matching/include/phat/compute_persistence_pairs.h | 11 | ||||
| -rw-r--r-- | matching/include/phat/helpers/misc.h | 3 |
7 files changed, 949 insertions, 45 deletions
diff --git a/matching/include/bifiltration.h b/matching/include/bifiltration.h index c9b64f1..2b4950f 100644 --- a/matching/include/bifiltration.h +++ b/matching/include/bifiltration.h @@ -38,7 +38,7 @@ namespace md { init(); } - DiagramKeeper weighted_slice_diagram(const DualPoint& line, int dim) const; + Diagram weighted_slice_diagram(const DualPoint& line, int dim) const; SimplexVector simplices() const { return simplices_; } @@ -100,8 +100,37 @@ namespace md { }; std::ostream& operator<<(std::ostream& os, const Bifiltration& bif); + + class BifiltrationProxy { + public: + BifiltrationProxy(const Bifiltration& bif, int dim = 0); + // return critical values of simplices that are important for current dimension (dim and dim+1) + PointVec positions() const; + // set current dimension + int set_dim(int new_dim); + + // wrappers of Bifiltration + int maximal_dim() const; + void translate(Real a); + Real minimal_coordinate() const; + Box bounding_box() const; + Real max_x() const; + Real max_y() const; + Real min_x() const; + Real min_y() const; + Diagram weighted_slice_diagram(const DualPoint& slice) const; + + private: + int dim_ { 0 }; + mutable PointVec cached_positions_; + Bifiltration bif_; + + void cache_positions() const; + }; } + + #endif //MATCHING_DISTANCE_BIFILTRATION_H //// The value type of OutputIterator is Simplex_in_2D_filtration diff --git a/matching/include/common_util.h b/matching/include/common_util.h index 674783b..b2ee22d 100644 --- a/matching/include/common_util.h +++ b/matching/include/common_util.h @@ -11,18 +11,51 @@ #include <map> #include "common_defs.h" +#include "phat/helpers/misc.h" namespace md { using Real = double; - using Index = long; + using RealVec = std::vector<Real>; + using Index = phat::index; + using IndexVec = std::vector<Index>; static constexpr Real pi = M_PI; using Column = std::vector<Index>; + struct IntPoint { + Index x; + Index y; + + IntPoint(Index x, Index y) :x(x), y(y) { } + + IntPoint() :x(0), y(0) { } + + inline bool operator==(const IntPoint& v) const + { + return this->x == v.x && this->y == v.y; + } + + // compare both coordinates, as needed in persistence + // do not overload operator<, requirements are not satisfied + inline bool is_less(const IntPoint& other, bool strict = false) const + { + if (x <= other.x and y <= other.y) { + if (strict) { + return x != other.x or y != other.y; + } + else { + return true; + } + } + return false; + } + }; + + struct Point { Real x; Real y; @@ -39,6 +72,12 @@ namespace md { y /= nor; } + inline void translate(Real a) + { + x += a; + y += a; + } + inline bool operator==(const Point& v) const { return this->x == v.x && this->y == v.y; @@ -60,6 +99,8 @@ namespace md { } }; + using PointVec = std::vector<Point>; + Point operator+(const Point& u, const Point& v); Point operator-(const Point& u, const Point& v); diff --git a/matching/include/matching_distance.h b/matching/include/matching_distance.h index 603b535..72743a4 100644 --- a/matching/include/matching_distance.h +++ b/matching/include/matching_distance.h @@ -6,6 +6,7 @@ #include <ostream> #include "spdlog/spdlog.h" +#include "spdlog/fmt/ostr.h" #include "common_defs.h" #include "cell_with_value.h" @@ -53,7 +54,7 @@ namespace md { struct CalculationParams { static constexpr int ALL_DIMENSIONS = -1; - Real hera_epsilon {0.01}; // relative error in hera call + Real hera_epsilon {0.001}; // relative error in hera call Real delta {0.1}; // relative error for matching distance int max_depth {6}; // maximal number of refinenemnts int initialization_depth {3}; @@ -76,9 +77,11 @@ namespace md { #endif }; + template<class DiagramProvider> class DistanceCalculator { - using DiagramProvider = md::Bifiltration; +// using DiagramProvider = md::Bifiltration; +// using DiagramProvider = md::ModulePresentation; using DualBox = md::DualBox; using CellValueVector = std::vector<CellWithValue>; @@ -90,29 +93,17 @@ namespace md { Real distance(); - int get_hera_calls_number() const - { - int result; - for(auto c : n_hera_calls_) result += c.second; - return result; - } - - int get_hera_calls_number(int dim) const; - - void clear_cache(); - -// for tests - make everything p + int get_hera_calls_number() const; +// for tests - make everything public // private: DiagramProvider module_a_; DiagramProvider module_b_; CalculationParams& params_; - int maximal_dim_ {0}; - - std::map<int, int> n_hera_calls_; + int n_hera_calls_; std::map<int, int> n_hera_calls_per_level_; - std::vector<Real> distances_; // indexed by dim + Real distance_; // if calculate_on_intermediate, then weighted distance // will be calculated on centers of each grid in between @@ -126,18 +117,18 @@ namespace md { Real get_max_y(int module) const; - void set_cell_central_value(CellWithValue& dual_cell, int dim); + void set_cell_central_value(CellWithValue& dual_cell); - Real distance_in_dimension(int dim); + Real get_distance(); - Real distance_in_dimension_pq(int dim); + Real get_distance_pq(); // temporary, to try priority queue Real get_max_possible_value(const CellWithValue* first_cell_ptr, int n_cells); - Real get_upper_bound(const CellWithValue& dual_cell, int dim, Real good_enough_upper_bound) const; + Real get_upper_bound(const CellWithValue& dual_cell, Real good_enough_upper_bound) const; - Real get_single_dgm_bound(const CellWithValue& dual_cell, ValuePoint vp, int module, int dim, + Real get_single_dgm_bound(const CellWithValue& dual_cell, ValuePoint vp, int module, Real good_enough_value) const; // this bound depends only on dual box @@ -155,14 +146,30 @@ namespace md { Real get_max_displacement_single_point(const CellWithValue& dual_cell, ValuePoint value_point, const Point& p) const; - void check_upper_bound(const CellWithValue& dual_cell, int dim) const; + void check_upper_bound(const CellWithValue& dual_cell) const; - Real distance_on_line(int dim, DualPoint line); - Real distance_on_line_const(int dim, DualPoint line) const; + Real distance_on_line(DualPoint line); + Real distance_on_line_const(DualPoint line) const; Real current_error(Real lower_bound, Real upper_bound); }; Real matching_distance(const Bifiltration& bif_a, const Bifiltration& bif_b, CalculationParams& params); + Real matching_distance(const ModulePresentation& mod_a, const ModulePresentation& mod_b, CalculationParams& params); + + // for upper bound experiment + struct UbExperimentRecord { + Real error; + Real lower_bound; + Real upper_bound; + CellWithValue cell; + long long int time; + long long int n_hera_calls; + }; + + std::ostream& operator<<(std::ostream& os, const UbExperimentRecord& r); + } + +#include "matching_distance.hpp" diff --git a/matching/include/matching_distance.hpp b/matching/include/matching_distance.hpp new file mode 100644 index 0000000..218f33b --- /dev/null +++ b/matching/include/matching_distance.hpp @@ -0,0 +1,784 @@ +namespace md { + + template<class K, class V> + void print_map(const std::map<K, V>& dic) + { + for(const auto kv : dic) { + fmt::print("{} -> {}\n", kv.first, kv.second); + } + } + + template<class T> + void DistanceCalculator<T>::check_upper_bound(const CellWithValue& dual_cell) const + { + spd::debug("Enter check_get_max_delta_on_cell"); + const int n_samples_lambda = 100; + const int n_samples_mu = 100; + DualBox db = dual_cell.dual_box(); + Real min_lambda = db.lambda_min(); + Real max_lambda = db.lambda_max(); + Real min_mu = db.mu_min(); + Real max_mu = db.mu_max(); + + Real h_lambda = (max_lambda - min_lambda) / n_samples_lambda; + Real h_mu = (max_mu - min_mu) / n_samples_mu; + for(int i = 1; i < n_samples_lambda; ++i) { + for(int j = 1; j < n_samples_mu; ++j) { + Real lambda = min_lambda + i * h_lambda; + Real mu = min_mu + j * h_mu; + DualPoint l(db.axis_type(), db.angle_type(), lambda, mu); + Real other_result = distance_on_line_const(l); + Real diff = fabs(dual_cell.stored_upper_bound() - other_result); + if (other_result > dual_cell.stored_upper_bound()) { + spd::error( + "in check_upper_bound, upper_bound = {}, other_result = {}, diff = {}\ndual_cell = {}", + dual_cell.stored_upper_bound(), other_result, diff, dual_cell); + throw std::runtime_error("Wrong delta estimate"); + } + } + } + spd::debug("Exit check_get_max_delta_on_cell"); + } + + // for all lines l, l' inside dual box, + // find the upper bound on the difference of weighted pushes of p + template<class T> + Real + DistanceCalculator<T>::get_max_displacement_single_point(const CellWithValue& dual_cell, ValuePoint vp, + const Point& p) const + { + assert(p.x >= 0 && p.y >= 0); + +#ifdef MD_DEBUG + std::vector<long long int> debug_ids = {3, 13, 54, 218, 350, 382, 484, 795, 2040, 8415, 44076}; + bool debug = false; // std::find(debug_ids.begin(), debug_ids.end(), dual_cell.id) != debug_ids.end(); +#endif + DualPoint line = dual_cell.value_point(vp); + const Real base_value = line.weighted_push(p); + + spd::debug("Enter get_max_displacement_single_point, p = {},\ndual_cell = {},\nline = {}, base_value = {}\n", p, + dual_cell, line, base_value); + + Real result = 0.0; + for(DualPoint dp : dual_cell.dual_box().critical_points(p)) { + Real dp_value = dp.weighted_push(p); + spd::debug( + "In get_max_displacement_single_point, p = {}, critical dp = {},\ndp_value = {}, diff = {},\ndual_cell = {}\n", + p, dp, dp_value, fabs(base_value - dp_value), dual_cell); + result = std::max(result, fabs(base_value - dp_value)); + } + +#ifdef MD_DO_FULL_CHECK + DualBox db = dual_cell.dual_box(); + std::uniform_real_distribution<Real> dlambda(db.lambda_min(), db.lambda_max()); + std::uniform_real_distribution<Real> dmu(db.mu_min(), db.mu_max()); + std::mt19937 gen(1); + for(int i = 0; i < 1000; ++i) { + Real lambda = dlambda(gen); + Real mu = dmu(gen); + DualPoint dp_random { db.axis_type(), db.angle_type(), lambda, mu }; + Real dp_value = dp_random.weighted_push(p); + if (fabs(base_value - dp_value) > result) { + spd::error("in get_max_displacement_single_point, p = {}, vp = {}\ndb = {}\nresult = {}, base_value = {}, dp_value = {}, dp_random = {}", + p, vp, db, result, base_value, dp_value, dp_random); + throw std::runtime_error("error in get_max_displacement_single_value"); + } + } +#endif + + return result; + } + + template<class T> + typename DistanceCalculator<T>::CellValueVector DistanceCalculator<T>::get_initial_dual_grid(Real& lower_bound) + { + CellValueVector result = get_refined_grid(params_.initialization_depth, false, true); + + lower_bound = -1.0; + for(const auto& dc : result) { + lower_bound = std::max(lower_bound, dc.max_corner_value()); + } + + assert(lower_bound >= 0); + + for(auto& dual_cell : result) { + Real good_enough_ub = get_good_enough_upper_bound(lower_bound); + Real max_value_on_cell = get_upper_bound(dual_cell, good_enough_ub); + dual_cell.set_max_possible_value(max_value_on_cell); + +#ifdef MD_DO_FULL_CHECK + check_upper_bound(dual_cell); +#endif + + spd::debug("DEBUG INIT: added cell {}", dual_cell); + } + + return result; + } + + template<class T> + typename DistanceCalculator<T>::CellValueVector + DistanceCalculator<T>::get_refined_grid(int init_depth, bool calculate_on_intermediate, bool calculate_on_last) + { + const Real y_max = std::max(module_a_.max_y(), module_b_.max_y()); + const Real x_max = std::max(module_a_.max_x(), module_b_.max_x()); + + const Real lambda_min = 0; + const Real lambda_max = 1; + + const Real mu_min = 0; + + DualBox x_flat(DualPoint(AxisType::x_type, AngleType::flat, lambda_min, mu_min), + DualPoint(AxisType::x_type, AngleType::flat, lambda_max, x_max)); + + DualBox x_steep(DualPoint(AxisType::x_type, AngleType::steep, lambda_min, mu_min), + DualPoint(AxisType::x_type, AngleType::steep, lambda_max, x_max)); + + DualBox y_flat(DualPoint(AxisType::y_type, AngleType::flat, lambda_min, mu_min), + DualPoint(AxisType::y_type, AngleType::flat, lambda_max, y_max)); + + DualBox y_steep(DualPoint(AxisType::y_type, AngleType::steep, lambda_min, mu_min), + DualPoint(AxisType::y_type, AngleType::steep, lambda_max, y_max)); + + CellWithValue x_flat_cell(x_flat, 0); + CellWithValue x_steep_cell(x_steep, 0); + CellWithValue y_flat_cell(y_flat, 0); + CellWithValue y_steep_cell(y_steep, 0); + + if (init_depth == 0) { + DualPoint diagonal_x_flat(AxisType::x_type, AngleType::flat, 1, 0); + + Real diagonal_value = distance_on_line(diagonal_x_flat); + n_hera_calls_per_level_[0]++; + + x_flat_cell.set_value_at(ValuePoint::lower_right, diagonal_value); + y_flat_cell.set_value_at(ValuePoint::lower_right, diagonal_value); + x_steep_cell.set_value_at(ValuePoint::lower_right, diagonal_value); + y_steep_cell.set_value_at(ValuePoint::lower_right, diagonal_value); + } + +#ifdef MD_DEBUG + x_flat_cell.id = 1; + x_steep_cell.id = 2; + y_flat_cell.id = 3; + y_steep_cell.id = 4; + CellWithValue::max_id = 4; +#endif + + CellValueVector result {x_flat_cell, x_steep_cell, y_flat_cell, y_steep_cell}; + + if (init_depth == 0) { + return result; + } + + CellValueVector refined_result; + + for(int i = 1; i <= init_depth; ++i) { + refined_result.clear(); + for(const auto& dual_cell : result) { + for(auto refined_cell : dual_cell.get_refined_cells()) { + // we calculate for init_dept - 1, not init_depth, + // because we want the cells to have value at a corner + if ((i == init_depth - 1 and calculate_on_last) or calculate_on_intermediate) + set_cell_central_value(refined_cell); + refined_result.push_back(refined_cell); + } + } + result = std::move(refined_result); + } + return result; + } + + template<class T> + DistanceCalculator<T>::DistanceCalculator(const T& a, + const T& b, + CalculationParams& params) + : + module_a_(a), + module_b_(b), + params_(params) + { + // make all coordinates non-negative + auto min_coord = std::min(module_a_.minimal_coordinate(), + module_b_.minimal_coordinate()); + if (min_coord < 0) { + module_a_.translate(-min_coord); + module_b_.translate(-min_coord); + } + + assert(std::min({module_a_.min_x(), module_b_.min_x(), module_a_.min_y(), + module_b_.min_y()}) >= 0); + + spd::info("DistanceCalculator constructed, module_a: max_x = {}, max_y = {}, module_b: max_x = {}, max_y = {}", + module_a_.max_x(), module_a_.max_y(), module_b_.max_x(), module_b_.max_y()); + } + + template<class T> + Real DistanceCalculator<T>::get_max_x(int module) const + { + return (module == 0) ? module_a_.max_x() : module_b_.max_x(); + } + + template<class T> + Real DistanceCalculator<T>::get_max_y(int module) const + { + return (module == 0) ? module_a_.max_y() : module_b_.max_y(); + } + + template<class T> + Real + DistanceCalculator<T>::get_local_refined_bound(const md::DualBox& dual_box) const + { + return get_local_refined_bound(0, dual_box) + get_local_refined_bound(1, dual_box); + } + + template<class T> + Real + DistanceCalculator<T>::get_local_refined_bound(int module, const md::DualBox& dual_box) const + { + spd::debug("Enter get_local_refined_bound, dual_box = {}", dual_box); + Real d_lambda = dual_box.lambda_max() - dual_box.lambda_min(); + Real d_mu = dual_box.mu_max() - dual_box.mu_min(); + Real result; + if (dual_box.axis_type() == AxisType::x_type) { + if (dual_box.is_flat()) { + result = dual_box.lambda_max() * d_mu + (get_max_x(module) - dual_box.mu_min()) * d_lambda; + } else { + result = d_mu + get_max_y(module) * d_lambda; + } + } else { + // y-type + if (dual_box.is_flat()) { + result = d_mu + get_max_x(module) * d_lambda; + } else { + // steep + result = dual_box.lambda_max() * d_mu + (get_max_y(module) - dual_box.mu_min()) * d_lambda; + } + } + return result; + } + + template<class T> + Real DistanceCalculator<T>::get_local_dual_bound(int module, const md::DualBox& dual_box) const + { + Real dlambda = dual_box.lambda_max() - dual_box.lambda_min(); + Real dmu = dual_box.mu_max() - dual_box.mu_min(); + + if (dual_box.is_flat()) { + return get_max_x(module) * dlambda + dmu; + } else { + return get_max_y(module) * dlambda + dmu; + } + } + + template<class T> + Real DistanceCalculator<T>::get_local_dual_bound(const md::DualBox& dual_box) const + { + return get_local_dual_bound(0, dual_box) + get_local_dual_bound(1, dual_box); + } + + template<class T> + Real DistanceCalculator<T>::get_upper_bound(const CellWithValue& dual_cell, Real good_enough_ub) const + { + assert(good_enough_ub >= 0); + + switch(params_.bound_strategy) { + case BoundStrategy::bruteforce: + return std::numeric_limits<Real>::max(); + + case BoundStrategy::local_dual_bound: + return dual_cell.min_value() + get_local_dual_bound(dual_cell.dual_box()); + + case BoundStrategy::local_dual_bound_refined: + return dual_cell.min_value() + get_local_refined_bound(dual_cell.dual_box()); + + case BoundStrategy::local_combined: { + Real cheap_upper_bound = dual_cell.min_value() + get_local_refined_bound(dual_cell.dual_box()); + if (cheap_upper_bound < good_enough_ub) { + return cheap_upper_bound; + } else { + [[fallthrough]]; + } + } + + case BoundStrategy::local_dual_bound_for_each_point: { + Real result = std::numeric_limits<Real>::max(); + for(ValuePoint vp : k_corner_vps) { + if (not dual_cell.has_value_at(vp)) { + continue; + } + + Real base_value = dual_cell.value_at(vp); + Real bound_dgm_a = get_single_dgm_bound(dual_cell, vp, 0, good_enough_ub); + + if (params_.stop_asap and bound_dgm_a + base_value >= good_enough_ub) { + // we want to return a valid upper bound, not just something that will prevent discarding the cell + // and we don't want to compute pushes for points in second bifiltration. + // so just return a constant time bound + return dual_cell.min_value() + get_local_refined_bound(dual_cell.dual_box()); + } + + Real bound_dgm_b = get_single_dgm_bound(dual_cell, vp, 1, + std::max(Real(0), good_enough_ub - bound_dgm_a)); + + result = std::min(result, base_value + bound_dgm_a + bound_dgm_b); + +#ifdef MD_DEBUG + spd::debug("In get_upper_bound, cell = {}", dual_cell); + spd::debug("In get_upper_bound, vp = {}, base_value = {}, bound_dgm_a = {}, bound_dgm_b = {}, result = {}", vp, base_value, bound_dgm_a, bound_dgm_b, result); +#endif + + if (params_.stop_asap and result < good_enough_ub) { + break; + } + } + return result; + } + } + // to suppress compiler warning + return std::numeric_limits<Real>::max(); + } + + // find maximal displacement of weighted points of m for all lines in dual_box + template<class T> + Real + DistanceCalculator<T>::get_single_dgm_bound(const CellWithValue& dual_cell, + ValuePoint vp, + int module, + [[maybe_unused]] Real good_enough_value) const + { + Real result = 0; + Point max_point; + + spd::debug( + "Enter get_single_dgm_bound, module = {}, dual_cell = {}, vp = {}, good_enough_value = {}, stop_asap = {}\n", + module, dual_cell, vp, good_enough_value, params_.stop_asap); + + const T& m = (module == 0) ? module_a_ : module_b_; + for(const auto& position : m.positions()) { + spd::debug("in get_single_dgm_bound, simplex = {}\n", position); + + Real x = get_max_displacement_single_point(dual_cell, vp, position); + + spd::debug("In get_single_dgm_bound, point = {}, displacement = {}", position, x); + + if (x > result) { + result = x; + max_point = position; + spd::debug("In get_single_dgm_bound, point = {}, result now = displacement = {}", position, x); + } + + if (params_.stop_asap and result > good_enough_value) { + // we want to return a valid upper bound, + // now we just see it is worse than we need, but it may be even more + // just return a valid upper bound + spd::debug("result {} > good_enough_value {}, exit and return refined bound {}", result, + good_enough_value, get_local_refined_bound(dual_cell.dual_box())); + result = get_local_refined_bound(dual_cell.dual_box()); + break; + } + } + + spd::debug("Exit get_single_dgm_bound,\ndual_cell = {}\nmodule = {}, result = {}, max_point = {}", dual_cell, + module, result, max_point); + + return result; + } + + template<class T> + Real DistanceCalculator<T>::distance() + { + return get_distance_pq(); + } + + // calculate weighted bottleneneck distance between slices on line + // increments hera calls counter + template<class T> + Real DistanceCalculator<T>::distance_on_line(DualPoint line) + { + ++n_hera_calls_; + Real result = distance_on_line_const(line); + return result; + } + + template<class T> + Real DistanceCalculator<T>::distance_on_line_const(DualPoint line) const + { + // TODO: think about this - how to call Hera + auto dgm_a = module_a_.weighted_slice_diagram(line); + auto dgm_b = module_b_.weighted_slice_diagram(line); + Real result; + if (params_.hera_epsilon > static_cast<Real>(0)) { + result = hera::bottleneckDistApprox(dgm_a, dgm_b, params_.hera_epsilon) / ( params_.hera_epsilon + 1); + } else { + result = hera::bottleneckDistExact(dgm_a, dgm_b); + } + if (n_hera_calls_ % 100 == 1) { + spd::debug("Calling Hera, dgm_a.size = {}, dgm_b.size = {}, line = {}, result = {}", dgm_a.size(), + dgm_b.size(), line, result); + } else { + spd::debug("Calling Hera, dgm_a.size = {}, dgm_b.size = {}, line = {}, result = {}", dgm_a.size(), + dgm_b.size(), line, result); + } + return result; + } + + template<class T> + Real DistanceCalculator<T>::get_good_enough_upper_bound(Real lower_bound) const + { + Real result; + // in upper_bound strategy we only prune cells if they cannot improve the lower bound, + // otherwise the experiment is supposed to run indefinitely + if (params_.traverse_strategy == TraverseStrategy::upper_bound) { + result = lower_bound; + } else { + result = (1.0 + params_.delta) * lower_bound; + } + return result; + } + + // helper function + // calculate weighted bt distance on cell center, + // assign distance value to cell, keep it in heat_map, and return + template<class T> + void DistanceCalculator<T>::set_cell_central_value(CellWithValue& dual_cell) + { + DualPoint central_line {dual_cell.center()}; + + spd::debug("In set_cell_central_value, processing dual cell = {}, line = {}", dual_cell.dual_box(), + central_line); + Real new_value = distance_on_line(central_line); + n_hera_calls_per_level_[dual_cell.level() + 1]++; + dual_cell.set_value_at(ValuePoint::center, new_value); + params_.actual_max_depth = std::max(params_.actual_max_depth, dual_cell.level() + 1); + +#ifdef PRINT_HEAT_MAP + if (params_.bound_strategy == BoundStrategy::bruteforce) { + spd::debug("In set_cell_central_value, adding to heat_map pair {} - {}", dual_cell.center(), new_value); + if (dual_cell.level() > params_.initialization_depth + 1 + and params_.heat_maps[dual_cell.level()].count(dual_cell.center()) > 0) { + auto existing = params_.heat_maps[dual_cell.level()].find(dual_cell.center()); + spd::debug("EXISTING: {} -> {}", existing->first, existing->second); + } + assert(dual_cell.level() <= params_.initialization_depth + 1 + or params_.heat_maps[dual_cell.level()].count(dual_cell.center()) == 0); + params_.heat_maps[dual_cell.level()][dual_cell.center()] = new_value; + } +#endif + } + + // quick-and-dirty hack to efficiently traverse priority queue with dual cells + // returns maximal possible value on all cells in queue + // assumes that the underlying container is vector! + // cell_ptr: pointer to the first element in queue + // n_cells: queue size + template<class T> + Real DistanceCalculator<T>::get_max_possible_value(const CellWithValue* cell_ptr, int n_cells) + { + Real result = (n_cells > 0) ? cell_ptr->stored_upper_bound() : 0; + for(int i = 0; i < n_cells; ++i, ++cell_ptr) { + result = std::max(result, cell_ptr->stored_upper_bound()); + } + return result; + } + + // helper function: + // return current error from lower and upper bounds + // and save it in params_ (hence not const) + template<class T> + Real DistanceCalculator<T>::current_error(Real lower_bound, Real upper_bound) + { + Real current_error = (lower_bound > 0.0) ? (upper_bound - lower_bound) / lower_bound + : std::numeric_limits<Real>::max(); + + params_.actual_error = current_error; + + if (current_error < params_.delta) { + spd::debug( + "Threshold achieved! bound_strategy = {}, traverse_strategy = {}, upper_bound = {}, current_error = {}", + params_.bound_strategy, params_.traverse_strategy, upper_bound, current_error); + } + return current_error; + } + + // return matching distance + // use priority queue to store dual cells + // comparison function depends on the strategies in params_ + // ressets hera calls counter + template<class T> + Real DistanceCalculator<T>::get_distance_pq() + { + std::map<int, long> n_cells_considered; + std::map<int, long> n_cells_pushed_into_queue; + long int n_too_deep_cells = 0; + std::map<int, long> n_cells_discarded; + std::map<int, long> n_cells_pruned; + + spd::info("Enter get_distance_pq, bound strategy = {}, traverse strategy = {}, stop_asap = {} ", + params_.bound_strategy, params_.traverse_strategy, params_.stop_asap); + + std::chrono::high_resolution_clock timer; + auto start_time = timer.now(); + + n_hera_calls_ = 0; + n_hera_calls_per_level_.clear(); + + + // if cell is too deep and is not pushed into queue, + // we still need to take its max value into account; + // the max over such cells is stored in max_result_on_too_fine_cells + Real upper_bound_on_deep_cells = -1; + + spd::debug("Started iterations in dual space, delta = {}, bound_strategy = {}", params_.delta, + params_.bound_strategy); + // user-defined less lambda function + // to regulate priority queue depending on strategy + auto dual_cell_less = [this](const CellWithValue& a, const CellWithValue& b) { + + int a_level = a.level(); + int b_level = b.level(); + Real a_value = a.max_corner_value(); + Real b_value = b.max_corner_value(); + Real a_ub = a.stored_upper_bound(); + Real b_ub = b.stored_upper_bound(); + if (this->params_.traverse_strategy == TraverseStrategy::upper_bound and + (not a.has_max_possible_value() or not b.has_max_possible_value())) { + throw std::runtime_error("no upper bound on cell"); + } + DualPoint a_lower_left = a.dual_box().lower_left(); + DualPoint b_lower_left = b.dual_box().lower_left(); + + switch(this->params_.traverse_strategy) { + // in both breadth_first searches we want coarser cells + // to be processed first. Cells with smaller level must be larger, + // hence the minus in front of level + case TraverseStrategy::breadth_first: + return std::make_tuple(-a_level, a_lower_left) + < std::make_tuple(-b_level, b_lower_left); + case TraverseStrategy::breadth_first_value: + return std::make_tuple(-a_level, a_value, a_lower_left) + < std::make_tuple(-b_level, b_value, b_lower_left); + case TraverseStrategy::depth_first: + return std::make_tuple(a_value, a_level, a_lower_left) + < std::make_tuple(b_value, b_level, b_lower_left); + case TraverseStrategy::upper_bound: + return std::make_tuple(a_ub, a_level, a_lower_left) + < std::make_tuple(b_ub, b_level, b_lower_left); + default: + throw std::runtime_error("Forgotten case"); + } + }; + + std::priority_queue<CellWithValue, CellValueVector, decltype(dual_cell_less)> dual_cells_queue( + dual_cell_less); + + // weighted bt distance on the center of current cell + Real lower_bound = std::numeric_limits<Real>::min(); + + // init pq and lower bound + for(auto& init_cell : get_initial_dual_grid(lower_bound)) { + dual_cells_queue.push(init_cell); + } + + Real upper_bound = get_max_possible_value(&dual_cells_queue.top(), dual_cells_queue.size()); + + std::vector<UbExperimentRecord> ub_experiment_results; + + while(not dual_cells_queue.empty()) { + + CellWithValue dual_cell = dual_cells_queue.top(); + dual_cells_queue.pop(); + assert(dual_cell.has_corner_value() + and dual_cell.has_max_possible_value() + and dual_cell.max_corner_value() <= upper_bound); + + n_cells_considered[dual_cell.level()]++; + + bool discard_cell = false; + + if (not params_.stop_asap) { + // if stop_asap is on, it is safer to never discard a cell + if (params_.bound_strategy == BoundStrategy::bruteforce) { + discard_cell = false; + } else if (params_.traverse_strategy == TraverseStrategy::upper_bound) { + discard_cell = (dual_cell.stored_upper_bound() <= lower_bound); + } else { + discard_cell = (dual_cell.stored_upper_bound() <= (1.0 + params_.delta) * lower_bound); + } + } + + spd::debug( + "CURRENT CELL bound_strategy = {}, traverse_strategy = {}, dual cell: {}, upper_bound = {}, lower_bound = {}, current_error = {}, discard_cell = {}", + params_.bound_strategy, params_.traverse_strategy, dual_cell, upper_bound, lower_bound, + current_error(lower_bound, upper_bound), discard_cell); + + if (discard_cell) { + n_cells_discarded[dual_cell.level()]++; + continue; + } + + // until now, dual_cell knows its value in one of its corners + // new_value will be the weighted distance at its center + set_cell_central_value(dual_cell); + Real new_value = dual_cell.value_at(ValuePoint::center); + lower_bound = std::max(new_value, lower_bound); + + spd::debug("Processed cell = {}, weighted value = {}, lower_bound = {}", dual_cell, new_value, lower_bound); + + assert(upper_bound >= lower_bound); + + if (current_error(lower_bound, upper_bound) < params_.delta) { + break; + } + + // refine cell and push 4 smaller cells into queue + for(auto refined_cell : dual_cell.get_refined_cells()) { + + if (refined_cell.num_values() == 0) + throw std::runtime_error("no value on cell"); + + // if delta is smaller than good_enough_value, it allows to prune cell + Real good_enough_ub = get_good_enough_upper_bound(lower_bound); + + // upper bound of the parent holds for refined_cell + // and can sometimes be smaller! + Real upper_bound_on_refined_cell = std::min(dual_cell.stored_upper_bound(), + get_upper_bound(refined_cell, good_enough_ub)); + + spd::debug("upper_bound_on_refined_cell = {}, dual_cell.stored_upper_bound = {}, get_upper_bound = {}", + upper_bound_on_refined_cell, dual_cell.stored_upper_bound(), + get_upper_bound(refined_cell, good_enough_ub)); + + refined_cell.set_max_possible_value(upper_bound_on_refined_cell); + +#ifdef MD_DO_FULL_CHECK + check_upper_bound(refined_cell); +#endif + + bool prune_cell = false; + + if (refined_cell.level() <= params_.max_depth) { + // cell might be added to queue; if it is not added, its maximal value can be safely ignored + if (params_.traverse_strategy == TraverseStrategy::upper_bound) { + prune_cell = (refined_cell.stored_upper_bound() <= lower_bound); + } else if (params_.bound_strategy != BoundStrategy::bruteforce) { + prune_cell = (refined_cell.stored_upper_bound() <= (1.0 + params_.delta) * lower_bound); + } + if (prune_cell) + n_cells_pruned[refined_cell.level()]++; +// prune_cell = (max_result_on_refined_cell <= lower_bound); + } else { + // cell is too deep, it won't be added to queue + // we must memorize maximal value on this cell, because we won't see it anymore + prune_cell = true; + if (refined_cell.stored_upper_bound() > (1 + params_.delta) * lower_bound) { + n_too_deep_cells++; + } + upper_bound_on_deep_cells = std::max(upper_bound_on_deep_cells, refined_cell.stored_upper_bound()); + } + + spd::debug( + "In get_distance_pq, loop over refined cells, bound_strategy = {}, traverse_strategy = {}, refined cell: {}, max_value_on_cell = {}, upper_bound = {}, current_error = {}, prune_cell = {}", + params_.bound_strategy, params_.traverse_strategy, refined_cell, + refined_cell.stored_upper_bound(), upper_bound, current_error(lower_bound, upper_bound), + prune_cell); + + if (not prune_cell) { + n_cells_pushed_into_queue[refined_cell.level()]++; + dual_cells_queue.push(refined_cell); + } + } // end loop over refined cells + + if (dual_cells_queue.empty()) + upper_bound = std::max(upper_bound, upper_bound_on_deep_cells); + else + upper_bound = std::max(upper_bound_on_deep_cells, + get_max_possible_value(&dual_cells_queue.top(), dual_cells_queue.size())); + + if (params_.traverse_strategy == TraverseStrategy::upper_bound) { + upper_bound = dual_cells_queue.top().stored_upper_bound(); + + if (get_hera_calls_number() < 20 || get_hera_calls_number() % 20 == 0) { + auto elapsed = timer.now() - start_time; + UbExperimentRecord ub_exp_record; + + ub_exp_record.error = current_error(lower_bound, upper_bound); + ub_exp_record.lower_bound = lower_bound; + ub_exp_record.upper_bound = upper_bound; + ub_exp_record.cell = dual_cells_queue.top(); + ub_exp_record.n_hera_calls = n_hera_calls_; + ub_exp_record.time = std::chrono::duration_cast<std::chrono::milliseconds>(elapsed).count(); + +#ifdef MD_DO_CHECKS + if (ub_experiment_results.size() > 0) { + auto prev = ub_experiment_results.back(); + if (upper_bound > prev.upper_bound) { + spd::error("ALARM 1, upper_bound = {}, top = {}, prev.ub = {}, prev cell = {}, lower_bound = {}, prev.lower_bound = {}", + upper_bound, ub_exp_record.cell, prev.upper_bound, prev.cell, lower_bound, prev.lower_bound); + throw std::runtime_error("die"); + } + + if (lower_bound < prev.lower_bound) { + spd::error("ALARM 2, lower_bound = {}, prev.lower_bound = {}, top = {}, prev.ub = {}, prev cell = {}", lower_bound, prev.lower_bound, ub_exp_record.cell, prev.upper_bound, prev.cell); + throw std::runtime_error("die"); + } + } +#endif + + ub_experiment_results.emplace_back(ub_exp_record); + + fmt::print(std::cerr, "[UB_EXPERIMENT]\t{}\n", ub_exp_record); + } + } + + assert(upper_bound >= lower_bound); + + if (current_error(lower_bound, upper_bound) < params_.delta) { + break; + } + } + + params_.actual_error = current_error(lower_bound, upper_bound); + + if (n_too_deep_cells > 0) { + spd::warn( + "Error not guaranteed, there were {} too deep cells. Actual error = {}. Increase max_depth or delta", + n_too_deep_cells, params_.actual_error); + } + // otherwise actual_error in params can be larger than delta, + // but this is OK + + spd::info("#############################################################"); + spd::info( + "Exiting get_distance_pq, bound_strategy = {}, traverse_strategy = {}, lower_bound = {}, upper_bound = {}, current_error = {}, actual_max_level = {}", + params_.bound_strategy, params_.traverse_strategy, lower_bound, + upper_bound, params_.actual_error, params_.actual_max_depth); + + spd::info("#############################################################"); + + bool print_stats = true; + if (print_stats) { + fmt::print("EXIT STATS, cells considered:\n"); + print_map(n_cells_considered); + fmt::print("EXIT STATS, cells discarded:\n"); + print_map(n_cells_discarded); + fmt::print("EXIT STATS, cells pruned:\n"); + print_map(n_cells_pruned); + fmt::print("EXIT STATS, cells pushed:\n"); + print_map(n_cells_pushed_into_queue); + fmt::print("EXIT STATS, hera calls:\n"); + print_map(n_hera_calls_per_level_); + + fmt::print("EXIT STATS, too deep cells with high value: {}\n", n_too_deep_cells); + } + + return lower_bound; + } + + template<class T> + int DistanceCalculator<T>::get_hera_calls_number() const + { + return n_hera_calls_; + } + +}
\ No newline at end of file diff --git a/matching/include/persistence_module.h b/matching/include/persistence_module.h index 50ae8d2..ea9cfd0 100644 --- a/matching/include/persistence_module.h +++ b/matching/include/persistence_module.h @@ -12,37 +12,68 @@ namespace md { - class Relation { + + class ModulePresentation { public: - Relation() {} - Relation(const Point& _pos) - :position_(_pos) { } + struct Relation { - Real get_x() const { return position_.x; } + Point position_; + IntPoint bigrade_; - Real get_y() const { return position_.y; } + Relation() {} + Relation(const Point& _pos) : position_(_pos) { } - private: - Point position_; - }; + Real get_x() const { return position_.x; } + Real get_y() const { return position_.y; } - class PersistenceModule { - public: - using Box = md::Box; + IndexVec components_; - PersistenceModule() { } + }; - PersistenceModule(const std::string& fname); // read from file + ModulePresentation() { } - Diagram slice_diagram(const DualPoint& line); + ModulePresentation(const std::string& fname); // read from file + Diagram weighted_slice_diagram(const DualPoint& line) const; + + // translate all points by vector (a,a) + void translate(Real a); + + // return minimal value of x- and y-coordinates + Real minimal_coordinate() const { return std::min(min_x(), min_x()); } + + // return box that contains all positions of all simplices Box bounding_box() const; + friend std::ostream& operator<<(std::ostream& os, const ModulePresentation& mp); + + Real max_x() const { return max_x_; } + + Real max_y() const { return max_y_; } + + Real min_x() const { return min_x_; } + + Real min_y() const { return min_y_; } + + PointVec positions() const; + private: - std::vector<Point> generators_; + + Real max_x_ { std::numeric_limits<Real>::max() }; + Real max_y_ { std::numeric_limits<Real>::max() }; + Real min_x_ { -std::numeric_limits<Real>::max() }; + Real min_y_ { -std::numeric_limits<Real>::max() }; + Box bounding_box_; + + PointVec generators_; std::vector<Relation> relations_; + PointVec positions_; + + void init_boundaries(); + void project_generators(const DualPoint& slice, IndexVec& sorted_indices, RealVec& projections) const; + void project_relations(const DualPoint& slice, IndexVec& sorted_indices, RealVec& projections) const; }; // template<typename R = double> diff --git a/matching/include/phat/compute_persistence_pairs.h b/matching/include/phat/compute_persistence_pairs.h index 48be65c..06f5372 100644 --- a/matching/include/phat/compute_persistence_pairs.h +++ b/matching/include/phat/compute_persistence_pairs.h @@ -97,14 +97,23 @@ namespace phat { ReductionAlgorithm reduce; reduce( boundary_matrix ); pairs.clear(); + std::set<index> max_indices; + // finite pairs for( index idx = 0; idx < boundary_matrix.get_num_cols(); idx++ ) { if( !boundary_matrix.is_empty( idx ) ) { index birth = boundary_matrix.get_max_index( idx ); + max_indices.insert(birth); index death = idx; pairs.append_pair( birth, death ); } } - } + // infinite pairs: column idx is 0, and row idx does not contain a lowest one + for( index idx = 0; idx < boundary_matrix.get_num_cols(); idx++ ) { + if(boundary_matrix.is_empty(idx) && max_indices.count(idx) == 0 ) { + pairs.append_pair( idx, k_infinity_index); + } + } + } template< typename ReductionAlgorithm, typename Representation > void compute_persistence_pairs_dualized( persistence_pairs& pairs, boundary_matrix< Representation >& boundary_matrix ) { diff --git a/matching/include/phat/helpers/misc.h b/matching/include/phat/helpers/misc.h index 3e1ed56..5a5c682 100644 --- a/matching/include/phat/helpers/misc.h +++ b/matching/include/phat/helpers/misc.h @@ -35,6 +35,7 @@ #include <cmath>
#include <cstdlib>
#include <iterator>
+#include <limits>
// VS2008 and below unfortunately do not support stdint.h
#if defined(_MSC_VER)&& _MSC_VER < 1600
@@ -55,6 +56,8 @@ namespace phat { typedef int64_t index;
typedef int8_t dimension;
typedef std::vector< index > column;
+
+ constexpr index k_infinity_index = std::numeric_limits<index>::max();
}
// OpenMP (proxy) functions
|