Skip to content

Commit

Permalink
Resolving conflict in PersistenceDiagramClustering
Browse files Browse the repository at this point in the history
  • Loading branch information
Keanu Sisouk authored and Keanu Sisouk committed Sep 22, 2023
1 parent baf3210 commit 51db9f0
Show file tree
Hide file tree
Showing 4 changed files with 214 additions and 210 deletions.
95 changes: 48 additions & 47 deletions core/base/persistenceDiagramClustering/PDBarycenter.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,14 +34,14 @@ void ttk::PDBarycenter::runMatching(
std::vector<std::vector<MatchingType>> *all_matchings,
bool use_kdt,
bool actual_distance) {
Timer time_matchings;
Timer const time_matchings;

double local_cost = *total_cost;
#ifdef TTK_ENABLE_OPENMP
#pragma omp parallel for num_threads(threadNumber_) schedule(dynamic, 1) reduction(+:local_cost)
#endif
for(int i = 0; i < numberOfInputs_; i++) {
double delta_lim = 0.01;
double const delta_lim = 0.01;
PersistenceDiagramAuction auction(
current_bidder_diagrams_[i], barycenter_goods_[i], wasserstein_,
geometrical_factor_, lambda_, delta_lim, kdt, correspondence_kdt_map,
Expand All @@ -56,15 +56,16 @@ void ttk::PDBarycenter::runMatching(
min_diag_price->at(i) = auction.getMinimalDiagonalPrice();
min_price->at(i) = getMinimalPrice(i);
std::vector<MatchingType> matchings;
double cost = auction.getMatchingsAndDistance(matchings, true);
double const cost = auction.getMatchingsAndDistance(matchings, true);
all_matchings->at(i) = matchings;
if(actual_distance) {
local_cost += sqrt(cost);
} else {
local_cost += cost;
}

double quotient = epsilon * auction.getAugmentedNumberOfBidders() / cost;
double const quotient
= epsilon * auction.getAugmentedNumberOfBidders() / cost;
precision_[i] = quotient < 1 ? 1. / sqrt(1 - quotient) - 1 : 10;
if(auction.getRelativePrecision() == 0)
precision_[i] = 0;
Expand Down Expand Up @@ -95,7 +96,7 @@ void ttk::PDBarycenter::runMatchingAuction(
geometrical_factor_, lambda_, 0.01, kdt, correspondence_kdt_map, 0,
(*min_diag_price)[i], use_kdt);
std::vector<MatchingType> matchings;
double cost = auction.run(matchings, i);
double const cost = auction.run(matchings, i);
all_matchings->at(i) = matchings;
if(actual_distance) {
local_cost += sqrt(cost);
Expand Down Expand Up @@ -140,7 +141,7 @@ std::vector<std::vector<ttk::MatchingType>> ttk::PDBarycenter::correctMatchings(
// 1. Invert the current_bidder_ids_ vector
std::vector<int> new_to_old_id(current_bidder_diagrams_[i].size());
for(size_t j = 0; j < current_bidder_ids_[i].size(); j++) {
int new_id = current_bidder_ids_[i][j];
int const new_id = current_bidder_ids_[i][j];
if(new_id >= 0) {
new_to_old_id[new_id] = j;
}
Expand All @@ -149,7 +150,7 @@ std::vector<std::vector<ttk::MatchingType>> ttk::PDBarycenter::correctMatchings(
std::vector<MatchingType> matchings_diagram_i;
for(size_t j = 0; j < previous_matchings[i].size(); j++) {
MatchingType m = previous_matchings[i][j];
int new_id = std::get<0>(m);
int const new_id = std::get<0>(m);
if(new_id >= 0 && std::get<1>(m) >= 0) {
std::get<0>(m) = new_to_old_id[new_id];
matchings_diagram_i.push_back(m);
Expand All @@ -163,10 +164,10 @@ std::vector<std::vector<ttk::MatchingType>> ttk::PDBarycenter::correctMatchings(
double ttk::PDBarycenter::updateBarycenter(
std::vector<std::vector<MatchingType>> &matchings) {
// 1. Initialize variables used in the sequel
Timer t_update;
size_t n_goods = barycenter_goods_[0].size();
Timer const t_update;
size_t const n_goods = barycenter_goods_[0].size();

size_t n_diagrams = current_bidder_diagrams_.size();
size_t const n_diagrams = current_bidder_diagrams_.size();
points_added_ = 0;
points_deleted_ = 0;
double max_shift = 0;
Expand Down Expand Up @@ -203,8 +204,8 @@ double ttk::PDBarycenter::updateBarycenter(
for(size_t j = 0; j < matchings.size(); j++) {
double weight = useCustomWeights_ ? (*customWeights_)[j] : 0.0;
for(size_t i = 0; i < matchings[j].size(); i++) {
int bidder_id = std::get<0>(matchings[j][i]);
int good_id = std::get<1>(matchings[j][i]);
int const bidder_id = std::get<0>(matchings[j][i]);
int const good_id = std::get<1>(matchings[j][i]);
if(good_id < 0 && bidder_id >= 0) {
// Future new barycenter point
points_to_append.emplace_back(
Expand Down Expand Up @@ -282,19 +283,19 @@ double ttk::PDBarycenter::updateBarycenter(
}

// TODO Weight by persistence
double new_crit_coord_x
double const new_crit_coord_x
= crit_coords_x[i] / (double)(n_diagrams - count_diag_matchings[i]);
double new_crit_coord_y
double const new_crit_coord_y
= crit_coords_y[i] / (double)(n_diagrams - count_diag_matchings[i]);
double new_crit_coord_z
double const new_crit_coord_z
= crit_coords_z[i] / (double)(n_diagrams - count_diag_matchings[i]);

// 3.3 Compute and store how much the point has shifted
// TODO adjust shift with geometrical_factor_
double dx = barycenter_goods_[0].at(i).x_ - new_x;
double dy = barycenter_goods_[0].at(i).y_ - new_y;
double shift = Geometry::pow(std::abs(dx), wasserstein_)
+ Geometry::pow(std::abs(dy), wasserstein_);
double const dx = barycenter_goods_[0].at(i).x_ - new_x;
double const dy = barycenter_goods_[0].at(i).y_ - new_y;
double const shift = Geometry::pow(std::abs(dx), wasserstein_)
+ Geometry::pow(std::abs(dy), wasserstein_);
if(shift > max_shift) {
max_shift = shift;
}
Expand Down Expand Up @@ -323,7 +324,7 @@ double ttk::PDBarycenter::updateBarycenter(
for(size_t i = 0; i < n_goods; i++) {
if(count_diag_matchings[i] == n_diagrams) {
points_deleted_ += 1;
double shift
double const shift
= 2
* Geometry::pow(
barycenter_goods_[0].at(i).getPersistence() / 2., wasserstein_);
Expand Down Expand Up @@ -360,7 +361,7 @@ double ttk::PDBarycenter::updateBarycenter(
std::get<2>(critical_coordinates));
}
barycenter_goods_[j].emplace_back(g);
double shift
double const shift
= 2
* Geometry::pow(
barycenter_goods_[j].at(g.id_).getPersistence() / 2., wasserstein_);
Expand Down Expand Up @@ -440,7 +441,7 @@ double ttk::PDBarycenter::enrichCurrentBidderDiagrams(
max_diagram_size = current_bidder_diagrams_[i].size();
}
}
int max_points_to_add = std::max(
int const max_points_to_add = std::max(
min_points_to_add, min_points_to_add + (int)(max_diagram_size / 10));

// 2. Get which points can be added, deduce the new minimal persistence
Expand All @@ -450,8 +451,8 @@ double ttk::PDBarycenter::enrichCurrentBidderDiagrams(

std::vector<double> persistences;
for(size_t j = 0; j < bidder_diagrams_[i].size(); j++) {
Bidder b = bidder_diagrams_[i].at(j);
double persistence = b.getPersistence();
Bidder const b = bidder_diagrams_[i].at(j);
double const persistence = b.getPersistence();
if(persistence >= min_persistence
&& persistence < previous_min_persistence) {
candidates_to_be_added[i].push_back(j);
Expand All @@ -463,9 +464,9 @@ double ttk::PDBarycenter::enrichCurrentBidderDiagrams(
return ((persistences[a] > persistences[b])
|| ((persistences[a] == persistences[b]) && (a > b)));
});
int size = candidates_to_be_added[i].size();
int const size = candidates_to_be_added[i].size();
if(size >= max_points_to_add) {
double last_persistence_added
double const last_persistence_added
= persistences[idx[i][max_points_to_add - 1]];
if(last_persistence_added > new_min_persistence) {
new_min_persistence = last_persistence_added;
Expand All @@ -479,7 +480,7 @@ double ttk::PDBarycenter::enrichCurrentBidderDiagrams(
int compteur_for_adding_points = 0;

for(int i = 0; i < numberOfInputs_; i++) {
int size = candidates_to_be_added[i].size();
int const size = candidates_to_be_added[i].size();
for(int j = 0; j < std::min(max_points_to_add, size); j++) {
Bidder b = bidder_diagrams_[i].at(candidates_to_be_added[i][idx[i][j]]);
if(b.getPersistence() >= new_min_persistence) {
Expand All @@ -491,7 +492,7 @@ double ttk::PDBarycenter::enrichCurrentBidderDiagrams(
current_bidder_ids_[i][candidates_to_be_added[i][idx[i][j]]]
= current_bidder_diagrams_[i].size() - 1;

int to_be_added_to_barycenter
int const to_be_added_to_barycenter
= deterministic_ ? compteur_for_adding_points % numberOfInputs_
: rand() % numberOfInputs_;
// We add the bidder as a good with probability 1/n_diagrams
Expand All @@ -516,8 +517,8 @@ double ttk::PDBarycenter::getMaxPersistence() {
BidderDiagram &D = bidder_diagrams_[i];
for(size_t j = 0; j < D.size(); j++) {
// Add bidder to bidders
Bidder &b = D.at(j);
double persistence = b.getPersistence();
Bidder const &b = D.at(j);
double const persistence = b.getPersistence();
if(persistence > max_persistence) {
max_persistence = persistence;
}
Expand All @@ -534,8 +535,8 @@ double ttk::PDBarycenter::getMinimalPrice(int i) {
return 0;
}
for(size_t j = 0; j < D.size(); j++) {
Good &b = D.at(j);
double price = b.getPrice();
Good const &b = D.at(j);
double const price = b.getPrice();
if(price < min_price) {
min_price = price;
}
Expand All @@ -552,8 +553,8 @@ double ttk::PDBarycenter::getLowestPersistence() {
BidderDiagram &D = bidder_diagrams_[i];
for(size_t j = 0; j < D.size(); j++) {
// Add bidder to bidders
Bidder &b = D.at(j);
double persistence = b.getPersistence();
Bidder const &b = D.at(j);
double const persistence = b.getPersistence();
if(persistence < lowest_persistence && persistence > 0) {
lowest_persistence = persistence;
}
Expand All @@ -579,7 +580,7 @@ void ttk::PDBarycenter::setInitialBarycenter(double min_persistence) {
int count = 0;
for(size_t j = 0; j < CTDiagram->size(); j++) {
// Add bidder to bidders
Good g = Good((*CTDiagram)[j], count, lambda_);
Good const g = Good((*CTDiagram)[j], count, lambda_);
if(g.getPersistence() >= min_persistence) {
goods.emplace_back(g);
count++;
Expand Down Expand Up @@ -617,7 +618,7 @@ typename ttk::PDBarycenter::KDTreePair ttk::PDBarycenter::getKDTree() const {
}

for(size_t idx = 0; idx < barycenter_goods_.size(); idx++) {
std::vector<double> empty_weights;
std::vector<double> const empty_weights;
weights.push_back(empty_weights);
for(size_t i = 0; i < barycenter_goods_[idx].size(); i++) {
const Good &g = barycenter_goods_[idx].at(i);
Expand All @@ -638,15 +639,15 @@ std::vector<std::vector<ttk::MatchingType>>
ttk::PDBarycenter::executeAuctionBarycenter(DiagramType &barycenter) {

std::vector<std::vector<MatchingType>> previous_matchings;
double min_persistence = 0;
double const min_persistence = 0;
double min_cost = std::numeric_limits<double>::max();
int last_min_cost_obtained = 0;

this->setBidderDiagrams();
this->setInitialBarycenter(
min_persistence); // false for a determinist initialization

double max_persistence = getMaxPersistence();
double const max_persistence = getMaxPersistence();

std::vector<double> min_diag_price(numberOfInputs_);
std::vector<double> min_price(numberOfInputs_);
Expand All @@ -655,7 +656,7 @@ std::vector<std::vector<ttk::MatchingType>>
min_price[i] = 0;
}

int min_points_to_add = std::numeric_limits<int>::max();
int const min_points_to_add = std::numeric_limits<int>::max();
this->enrichCurrentBidderDiagrams(2 * max_persistence, min_persistence,
min_diag_price, min_price,
min_points_to_add, false);
Expand All @@ -665,7 +666,7 @@ std::vector<std::vector<ttk::MatchingType>>
double total_cost;

while(!finished) {
Timer tm;
Timer const tm;

std::pair<std::unique_ptr<KDT>, std::vector<KDT *>> pair;
bool use_kdt = false;
Expand All @@ -686,13 +687,13 @@ std::vector<std::vector<ttk::MatchingType>>

barycenter.clear();
for(size_t j = 0; j < barycenter_goods_[0].size(); j++) {
Good &g = barycenter_goods_[0].at(j);
Good const &g = barycenter_goods_[0].at(j);
barycenter.emplace_back(PersistencePair{CriticalVertex{0, nt1_, g.x_, {}},
CriticalVertex{0, nt2_, g.y_, {}},
diagramType_, true});
}

bool actual_distance = (numberOfInputs_ == 2);
bool const actual_distance = (numberOfInputs_ == 2);
runMatchingAuction(&total_cost, sizes, *pair.first, pair.second,
&min_diag_price, &all_matchings, use_kdt,
actual_distance);
Expand Down Expand Up @@ -737,7 +738,7 @@ std::vector<std::vector<ttk::MatchingType>>
}
barycenter.resize(0);
for(size_t j = 0; j < barycenter_goods_[0].size(); j++) {
Good &g = barycenter_goods_[0].at(j);
Good const &g = barycenter_goods_[0].at(j);
barycenter.emplace_back(PersistencePair{CriticalVertex{0, nt1_, g.x_, {}},
CriticalVertex{0, nt2_, g.y_, {}},
diagramType_, true});
Expand All @@ -755,10 +756,10 @@ double ttk::PDBarycenter::computeRealCost() {
for(int i = 0; i < numberOfInputs_; i++) {
PersistenceDiagramAuction auction(
wasserstein_, geometrical_factor_, lambda_, 0.01, true);
GoodDiagram current_barycenter = barycenter_goods_[0];
BidderDiagram current_bidder_diagram = bidder_diagrams_[i];
GoodDiagram const current_barycenter = barycenter_goods_[0];
BidderDiagram const current_bidder_diagram = bidder_diagrams_[i];
auction.BuildAuctionDiagrams(current_bidder_diagram, current_barycenter);
double cost = auction.run(fake_matchings);
double const cost = auction.run(fake_matchings);
total_real_cost += cost * cost;
}
return sqrt(total_real_cost);
Expand All @@ -773,7 +774,7 @@ bool ttk::PDBarycenter::isPrecisionObjectiveMet(double precision_objective,
}
}
} else if(mode == 1) { // AVERAGE PRECISION
double average_precision
double const average_precision
= std::accumulate(precision_.begin(), precision_.end(), 0.0)
/ numberOfInputs_;
if(average_precision > precision_objective) {
Expand Down
Loading

0 comments on commit 51db9f0

Please sign in to comment.