TTK
Loading...
Searching...
No Matches
PDBarycenter.cpp
Go to the documentation of this file.
1
14
16
17#include <cmath>
18#include <cstdlib> /* srand, rand */
19#include <numeric>
20
21std::vector<std::vector<ttk::MatchingType>>
23 return executeAuctionBarycenter(barycenter);
24}
25
27 double *total_cost,
28 double epsilon,
29 std::vector<int> &sizes,
30 KDT &kdt,
31 std::vector<KDT *> &correspondence_kdt_map,
32 std::vector<double> *min_diag_price,
33 std::vector<double> *min_price,
34 std::vector<std::vector<MatchingType>> *all_matchings,
35 bool use_kdt,
36 bool actual_distance) {
37 Timer time_matchings;
38
39 double local_cost = *total_cost;
40#ifdef TTK_ENABLE_OPENMP
41#pragma omp parallel for num_threads(threadNumber_) schedule(dynamic, 1) reduction(+:local_cost)
42#endif
43 for(int i = 0; i < numberOfInputs_; i++) {
44 double delta_lim = 0.01;
46 current_bidder_diagrams_[i], barycenter_goods_[i], wasserstein_,
47 geometrical_factor_, lambda_, delta_lim, kdt, correspondence_kdt_map,
48 epsilon, min_diag_price->at(i), use_kdt);
49 int n_biddings = 0;
50 auction.initLowerBoundCostWeight(delta_lim);
51 auction.initLowerBoundCost(i);
52 auction.buildUnassignedBidders();
53 auction.reinitializeGoods();
54 auction.runAuctionRound(n_biddings, i);
55 auction.updateDiagonalPrices();
56 min_diag_price->at(i) = auction.getMinimalDiagonalPrice();
57 min_price->at(i) = getMinimalPrice(i);
58 std::vector<MatchingType> matchings;
59 double cost = auction.getMatchingsAndDistance(matchings, true);
60 all_matchings->at(i) = matchings;
61 if(actual_distance) {
62 local_cost += sqrt(cost);
63 } else {
64 local_cost += cost;
65 }
66
67 double quotient = epsilon * auction.getAugmentedNumberOfBidders() / cost;
68 precision_[i] = quotient < 1 ? 1. / sqrt(1 - quotient) - 1 : 10;
69 if(auction.getRelativePrecision() == 0)
70 precision_[i] = 0;
71 // Resizes the diagram which was enrich with diagonal bidders
72 // during the auction
73 // TODO do this inside the auction !
74 current_bidder_diagrams_[i].resize(sizes[i]);
75 }
76 *total_cost = local_cost;
77}
78
80 double *total_cost,
81 std::vector<int> &sizes,
82 KDT &kdt,
83 std::vector<KDT *> &correspondence_kdt_map,
84 std::vector<double> *min_diag_price,
85 std::vector<std::vector<MatchingType>> *all_matchings,
86 bool use_kdt,
87 bool actual_distance) {
88 double local_cost = *total_cost;
89#ifdef TTK_ENABLE_OPENMP
90#pragma omp parallel for num_threads(threadNumber_) schedule(dynamic, 1) reduction(+:local_cost)
91#endif
92 for(int i = 0; i < numberOfInputs_; i++) {
94 current_bidder_diagrams_[i], barycenter_goods_[i], wasserstein_,
95 geometrical_factor_, lambda_, 0.01, kdt, correspondence_kdt_map, 0,
96 (*min_diag_price)[i], use_kdt);
97 std::vector<MatchingType> matchings;
98 double cost = auction.run(matchings, i);
99 all_matchings->at(i) = matchings;
100 if(actual_distance) {
101 local_cost += sqrt(cost);
102 } else {
103 local_cost += cost;
104 }
105
106 // Resizes the diagram which was enrich with diagonal bidders
107 // during the auction
108 // TODO do this inside the auction !
109 current_bidder_diagrams_[i].resize(sizes[i]);
110 }
111 *total_cost = local_cost;
112}
113
115 std::vector<std::vector<MatchingType>> &matchings,
116 std::vector<std::vector<MatchingType>> &previous_matchings) {
117 if(points_added_ > 0 || points_deleted_ > 0 || previous_matchings.empty()) {
118 return false;
119 }
120
121 for(size_t j = 0; j < matchings.size(); j++) {
122 for(size_t i = 0; i < matchings[j].size(); i++) {
123 MatchingType t = matchings[j][i];
124 MatchingType previous_t = previous_matchings[j][i];
125
126 if(std::get<1>(t) != std::get<1>(previous_t)
127 && (std::get<0>(t) >= 0 && std::get<0>(previous_t) >= 0)) {
128 return false;
129 }
130 }
131 }
132 return true;
133}
134
135std::vector<std::vector<ttk::MatchingType>> ttk::PDBarycenter::correctMatchings(
136 std::vector<std::vector<MatchingType>> &previous_matchings) {
137
138 std::vector<std::vector<MatchingType>> corrected_matchings(numberOfInputs_);
139 for(int i = 0; i < numberOfInputs_; i++) {
140 // 1. Invert the current_bidder_ids_ vector
141 std::vector<int> new_to_old_id(current_bidder_diagrams_[i].size());
142 for(size_t j = 0; j < current_bidder_ids_[i].size(); j++) {
143 int new_id = current_bidder_ids_[i][j];
144 if(new_id >= 0) {
145 new_to_old_id[new_id] = j;
146 }
147 }
148 // 2. Reconstruct the matchings
149 std::vector<MatchingType> matchings_diagram_i;
150 for(size_t j = 0; j < previous_matchings[i].size(); j++) {
151 MatchingType m = previous_matchings[i][j];
152 int new_id = std::get<0>(m);
153 if(new_id >= 0 && std::get<1>(m) >= 0) {
154 std::get<0>(m) = new_to_old_id[new_id];
155 matchings_diagram_i.push_back(m);
156 }
157 }
158 corrected_matchings[i] = matchings_diagram_i;
159 }
160 return corrected_matchings;
161}
162
164 std::vector<std::vector<MatchingType>> &matchings) {
165 // 1. Initialize variables used in the sequel
166 Timer t_update;
167 size_t n_goods = barycenter_goods_[0].size();
168
169 size_t n_diagrams = current_bidder_diagrams_.size();
170 points_added_ = 0;
171 points_deleted_ = 0;
172 double max_shift = 0;
173
174 std::vector<size_t> count_diag_matchings(
175 n_goods); // Number of diagonal matchings for each point of the barycenter
176 std::vector<double> x(n_goods);
177 std::vector<double> y(n_goods);
178 std::vector<double> crit_coords_x(n_goods);
179 std::vector<double> crit_coords_y(n_goods);
180 std::vector<double> crit_coords_z(n_goods);
181 for(size_t i = 0; i < n_goods; i++) {
182 count_diag_matchings[i] = 0;
183 x[i] = 0;
184 y[i] = 0;
185 crit_coords_x[i] = 0;
186 crit_coords_y[i] = 0;
187 crit_coords_z[i] = 0;
188 }
189 std::vector<double> min_prices(n_diagrams);
190 for(size_t j = 0; j < n_diagrams; j++) {
191 min_prices[j] = std::numeric_limits<double>::max();
192 }
193
194 std::vector<Bidder *>
195 points_to_append; // Will collect bidders linked to diagonal
196 // 2. Preprocess the matchings
197 for(size_t j = 0; j < matchings.size(); j++) {
198 for(size_t i = 0; i < matchings[j].size(); i++) {
199 int bidder_id = std::get<0>(matchings[j][i]);
200 int good_id = std::get<1>(matchings[j][i]);
201 if(good_id < 0 && bidder_id >= 0) {
202 // Future new barycenter point
203 points_to_append.push_back(&current_bidder_diagrams_[j].at(bidder_id));
204 }
205
206 else if(good_id >= 0 && bidder_id >= 0) {
207 // Update coordinates (to be divided by the number of diagrams later on)
208 x[good_id] += current_bidder_diagrams_[j].at(bidder_id).x_;
209 y[good_id] += current_bidder_diagrams_[j].at(bidder_id).y_;
210 if(geometrical_factor_ < 1) {
211 const auto &critical_coordinates = current_bidder_diagrams_[j]
212 .at(bidder_id)
213 .GetCriticalCoordinates();
214 crit_coords_x[good_id] += critical_coordinates[0];
215 crit_coords_y[good_id] += critical_coordinates[1];
216 crit_coords_z[good_id] += critical_coordinates[2];
217 }
218 } else if(good_id >= 0 && bidder_id < 0) {
219 // Counting the number of times this barycenter point is linked to the
220 // diagonal
221 count_diag_matchings[good_id] = count_diag_matchings[good_id] + 1;
222 }
223 }
224 }
225
226 // 3. Update the previous points of the barycenter
227 for(size_t i = 0; i < n_goods; i++) {
228 if(count_diag_matchings[i] < n_diagrams) {
229 // Barycenter point i is matched at least to one off-diagonal bidder
230 // 3.1 Compute the arithmetic mean of off-diagonal bidders linked to it
231 double x_bar = x[i] / (double)(n_diagrams - count_diag_matchings[i]);
232 double y_bar = y[i] / (double)(n_diagrams - count_diag_matchings[i]);
233 // 3.2 Compute the new coordinates of the point (the more linked to the
234 // diagonal it was, the closer to the diagonal it'll be)
235 double new_x = ((double)(n_diagrams - count_diag_matchings[i]) * x_bar
236 + (double)count_diag_matchings[i] * (x_bar + y_bar) / 2.)
237 / (double)n_diagrams;
238 double new_y = ((double)(n_diagrams - count_diag_matchings[i]) * y_bar
239 + (double)count_diag_matchings[i] * (x_bar + y_bar) / 2.)
240 / (double)n_diagrams;
241 // TODO Weight by persistence
242 double new_crit_coord_x
243 = crit_coords_x[i] / (double)(n_diagrams - count_diag_matchings[i]);
244 double new_crit_coord_y
245 = crit_coords_y[i] / (double)(n_diagrams - count_diag_matchings[i]);
246 double new_crit_coord_z
247 = crit_coords_z[i] / (double)(n_diagrams - count_diag_matchings[i]);
248
249 // 3.3 Compute and store how much the point has shifted
250 // TODO adjust shift with geometrical_factor_
251 double dx = barycenter_goods_[0].at(i).x_ - new_x;
252 double dy = barycenter_goods_[0].at(i).y_ - new_y;
253 double shift = Geometry::pow(std::abs(dx), wasserstein_)
254 + Geometry::pow(std::abs(dy), wasserstein_);
255 if(shift > max_shift) {
256 max_shift = shift;
257 }
258 // 3.4 Update the position of the point
259 for(size_t j = 0; j < n_diagrams; j++) {
260 barycenter_goods_[j].at(i).SetCoordinates(new_x, new_y);
261 if(geometrical_factor_ < 1) {
262 barycenter_goods_[j].at(i).SetCriticalCoordinates(
263 new_crit_coord_x, new_crit_coord_y, new_crit_coord_z);
264 }
265 if(barycenter_goods_[j].at(i).getPrice() < min_prices[j]) {
266 min_prices[j] = barycenter_goods_[j].at(i).getPrice();
267 }
268 }
269 // TODO Reinitialize/play with prices here if you wish
270 }
271 }
272 for(size_t j = 0; j < n_diagrams; j++) {
273 if(min_prices[j] >= std::numeric_limits<double>::max() / 2.) {
274 min_prices[j] = 0;
275 }
276 }
277
278 // 4. Delete off-diagonal barycenter points not linked to any
279 // off-diagonal bidder
280 for(size_t i = 0; i < n_goods; i++) {
281 if(count_diag_matchings[i] == n_diagrams) {
282 points_deleted_ += 1;
283 double shift
284 = 2
286 barycenter_goods_[0].at(i).getPersistence() / 2., wasserstein_);
287 if(shift > max_shift) {
288 max_shift = shift;
289 }
290 for(size_t j = 0; j < n_diagrams; j++) {
291 barycenter_goods_[j].at(i).id_ = -1;
292 }
293 }
294 }
295
296 // 5. Append the new points to the barycenter
297 for(size_t k = 0; k < points_to_append.size(); k++) {
298 points_added_ += 1;
299 Bidder *b = points_to_append[k];
300 double gx
301 = (b->x_ + (n_diagrams - 1) * (b->x_ + b->y_) / 2.) / (n_diagrams);
302 double gy
303 = (b->y_ + (n_diagrams - 1) * (b->x_ + b->y_) / 2.) / (n_diagrams);
304 const auto &critical_coordinates = b->GetCriticalCoordinates();
305 for(size_t j = 0; j < n_diagrams; j++) {
306 Good g = Good(gx, gy, false, barycenter_goods_[j].size());
307 g.setPrice(min_prices[j]);
308 if(geometrical_factor_ < 1) {
309 g.SetCriticalCoordinates(std::get<0>(critical_coordinates),
310 std::get<1>(critical_coordinates),
311 std::get<2>(critical_coordinates));
312 }
313 barycenter_goods_[j].emplace_back(g);
314 double shift
315 = 2
317 barycenter_goods_[j].at(g.id_).getPersistence() / 2., wasserstein_);
318 if(shift > max_shift) {
319 max_shift = shift;
320 }
321 }
322 }
323
324 // 6. Finally, recreate barycenter_goods
325 for(size_t j = 0; j < n_diagrams; j++) {
326 int count = 0;
327 GoodDiagram new_barycenter;
328 for(size_t i = 0; i < barycenter_goods_[j].size(); i++) {
329 Good g = barycenter_goods_[j].at(i);
330 if(g.id_ != -1) {
331 g.id_ = count;
332 new_barycenter.emplace_back(g);
333 count++;
334 }
335 }
336 barycenter_goods_[j] = new_barycenter;
337 }
338
339 return max_shift;
340}
341
343 return rho * rho / 8.0;
344}
345
346double ttk::PDBarycenter::getRho(double epsilon) {
347 return std::sqrt(8.0 * epsilon);
348}
349
351
352 for(int i = 0; i < numberOfInputs_; i++) {
353 DiagramType *CTDiagram = &((*inputDiagrams_)[i]);
354
355 BidderDiagram bidders;
356 for(size_t j = 0; j < CTDiagram->size(); j++) {
357 // Add bidder to bidders
358 Bidder b((*CTDiagram)[j], j, lambda_);
359
360 b.setPositionInAuction(bidders.size());
361 bidders.emplace_back(b);
362 if(b.isDiagonal() || b.x_ == b.y_) {
363 this->printWrn("Diagonal point in diagram !!!");
364 }
365 }
366 bidder_diagrams_.push_back(bidders);
367 current_bidder_diagrams_.emplace_back();
368 std::vector<int> ids(bidders.size());
369 for(size_t j = 0; j < ids.size(); j++) {
370 ids[j] = -1;
371 }
372 current_bidder_ids_.push_back(ids);
373 }
374}
375
377 double previous_min_persistence,
378 double min_persistence,
379 std::vector<double> &initial_diagonal_prices,
380 std::vector<double> &initial_off_diagonal_prices,
381 int min_points_to_add,
382 bool add_points_to_barycenter) {
383
384 double new_min_persistence = min_persistence;
385
386 // 1. Get size of the largest current diagram, deduce the maximal number of
387 // points to append
388 size_t max_diagram_size{};
389 for(int i = 0; i < numberOfInputs_; i++) {
390 if(current_bidder_diagrams_[i].size() > max_diagram_size) {
391 max_diagram_size = current_bidder_diagrams_[i].size();
392 }
393 }
394 int max_points_to_add = std::max(
395 min_points_to_add, min_points_to_add + (int)(max_diagram_size / 10));
396
397 // 2. Get which points can be added, deduce the new minimal persistence
398 std::vector<std::vector<int>> candidates_to_be_added(numberOfInputs_);
399 std::vector<std::vector<int>> idx(numberOfInputs_);
400 for(int i = 0; i < numberOfInputs_; i++) {
401
402 std::vector<double> persistences;
403 for(size_t j = 0; j < bidder_diagrams_[i].size(); j++) {
404 Bidder b = bidder_diagrams_[i].at(j);
405 double persistence = b.getPersistence();
406 if(persistence >= min_persistence
407 && persistence < previous_min_persistence) {
408 candidates_to_be_added[i].push_back(j);
409 idx[i].push_back(idx[i].size());
410 persistences.push_back(persistence);
411 }
412 }
413 sort(idx[i].begin(), idx[i].end(), [&persistences](int &a, int &b) {
414 return ((persistences[a] > persistences[b])
415 || ((persistences[a] == persistences[b]) && (a > b)));
416 });
417 int size = candidates_to_be_added[i].size();
418 if(size >= max_points_to_add) {
419 double last_persistence_added
420 = persistences[idx[i][max_points_to_add - 1]];
421 if(last_persistence_added > new_min_persistence) {
422 new_min_persistence = last_persistence_added;
423 }
424 }
425 }
426
427 // 3. Add the points to the current diagrams
428
429 // only to give determinism
430 int compteur_for_adding_points = 0;
431
432 for(int i = 0; i < numberOfInputs_; i++) {
433 int size = candidates_to_be_added[i].size();
434 for(int j = 0; j < std::min(max_points_to_add, size); j++) {
435 Bidder b = bidder_diagrams_[i].at(candidates_to_be_added[i][idx[i][j]]);
436 if(b.getPersistence() >= new_min_persistence) {
437 b.id_ = current_bidder_diagrams_[i].size();
438 b.setPositionInAuction(current_bidder_diagrams_[i].size());
439 b.setDiagonalPrice(initial_diagonal_prices[i]);
440 current_bidder_diagrams_[i].emplace_back(b);
441 // b.id_ --> position of b in current_bidder_diagrams_[i]
442 current_bidder_ids_[i][candidates_to_be_added[i][idx[i][j]]]
443 = current_bidder_diagrams_[i].size() - 1;
444
445 int to_be_added_to_barycenter
446 = deterministic_ ? compteur_for_adding_points % numberOfInputs_
447 : rand() % numberOfInputs_;
448 // We add the bidder as a good with probability 1/n_diagrams
449 if(to_be_added_to_barycenter == 0 && add_points_to_barycenter) {
450 for(int k = 0; k < numberOfInputs_; k++) {
451 Good g = Good(b.x_, b.y_, false, barycenter_goods_[k].size());
452 g.setPrice(initial_off_diagonal_prices[k]);
453 g.SetCriticalCoordinates(b.coords_[0], b.coords_[1], b.coords_[2]);
454 barycenter_goods_[k].emplace_back(g);
455 }
456 }
457 }
458 compteur_for_adding_points++;
459 }
460 }
461 return new_min_persistence;
462}
463
465 double max_persistence = 0;
466 for(int i = 0; i < numberOfInputs_; i++) {
467 BidderDiagram &D = bidder_diagrams_[i];
468 for(size_t j = 0; j < D.size(); j++) {
469 // Add bidder to bidders
470 Bidder &b = D.at(j);
471 double persistence = b.getPersistence();
472 if(persistence > max_persistence) {
473 max_persistence = persistence;
474 }
475 }
476 }
477 return max_persistence;
478}
479
481 double min_price = std::numeric_limits<double>::max();
482
483 GoodDiagram &D = barycenter_goods_[i];
484 if(D.empty()) {
485 return 0;
486 }
487 for(size_t j = 0; j < D.size(); j++) {
488 Good &b = D.at(j);
489 double price = b.getPrice();
490 if(price < min_price) {
491 min_price = price;
492 }
493 }
494 if(min_price >= std::numeric_limits<double>::max() / 2.) {
495 return 0;
496 }
497 return min_price;
498}
499
501 double lowest_persistence = std::numeric_limits<double>::max();
502 for(int i = 0; i < numberOfInputs_; i++) {
503 BidderDiagram &D = bidder_diagrams_[i];
504 for(size_t j = 0; j < D.size(); j++) {
505 // Add bidder to bidders
506 Bidder &b = D.at(j);
507 double persistence = b.getPersistence();
508 if(persistence < lowest_persistence && persistence > 0) {
509 lowest_persistence = persistence;
510 }
511 }
512 }
513 if(lowest_persistence >= std::numeric_limits<double>::max() / 2.) {
514 return 0;
515 }
516 return lowest_persistence;
517}
518
519void ttk::PDBarycenter::setInitialBarycenter(double min_persistence) {
520 int size = 0;
521 int random_idx;
522 DiagramType *CTDiagram;
523 int iter = 0;
524 while(size == 0) {
525 random_idx
526 = deterministic_ ? iter % numberOfInputs_ : rand() % numberOfInputs_;
527 CTDiagram = &((*inputDiagrams_)[random_idx]);
528 for(int i = 0; i < numberOfInputs_; i++) {
529 GoodDiagram goods;
530 int count = 0;
531 for(size_t j = 0; j < CTDiagram->size(); j++) {
532 // Add bidder to bidders
533 Good g = Good((*CTDiagram)[j], count, lambda_);
534 if(g.getPersistence() >= min_persistence) {
535 goods.emplace_back(g);
536 count++;
537 }
538 }
539 if(static_cast<int>(barycenter_goods_.size()) < (i + 1)) {
540 barycenter_goods_.push_back(goods);
541 } else {
542 barycenter_goods_[i] = goods;
543 }
544 }
545 size = barycenter_goods_[0].size();
546 iter++;
547 }
548}
549
551 Timer tm;
552 auto kdt = std::make_unique<KDT>(true, wasserstein_);
553
554 const int dimension = geometrical_factor_ >= 1 ? 2 : 5;
555
556 std::vector<double> coordinates;
557 std::vector<std::vector<double>> weights;
558
559 for(size_t i = 0; i < barycenter_goods_[0].size(); i++) {
560 const Good &g = barycenter_goods_[0].at(i);
561 coordinates.push_back(geometrical_factor_ * g.x_);
562 coordinates.push_back(geometrical_factor_ * g.y_);
563 if(geometrical_factor_ < 1) {
564 coordinates.push_back((1 - geometrical_factor_) * g.coords_[0]);
565 coordinates.push_back((1 - geometrical_factor_) * g.coords_[1]);
566 coordinates.push_back((1 - geometrical_factor_) * g.coords_[2]);
567 }
568 }
569
570 for(size_t idx = 0; idx < barycenter_goods_.size(); idx++) {
571 std::vector<double> empty_weights;
572 weights.push_back(empty_weights);
573 for(size_t i = 0; i < barycenter_goods_[idx].size(); i++) {
574 const Good &g = barycenter_goods_[idx].at(i);
575 weights[idx].push_back(g.getPrice());
576 }
577 }
578 // Correspondence map : position in barycenter_goods_ --> KDT node
579
580 auto correspondence_kdt_map
581 = kdt->build(coordinates.data(), barycenter_goods_[0].size(), dimension,
582 weights, barycenter_goods_.size());
583 this->printMsg(" Building KDTree", 1, tm.getElapsedTime(),
585 return std::make_pair(std::move(kdt), correspondence_kdt_map);
586}
587
588std::vector<std::vector<ttk::MatchingType>>
590
591 std::vector<std::vector<MatchingType>> previous_matchings;
592 double min_persistence = 0;
593 double min_cost = std::numeric_limits<double>::max();
594 int last_min_cost_obtained = 0;
595
596 this->setBidderDiagrams();
597 this->setInitialBarycenter(
598 min_persistence); // false for a determinist initialization
599
600 double max_persistence = getMaxPersistence();
601
602 std::vector<double> min_diag_price(numberOfInputs_);
603 std::vector<double> min_price(numberOfInputs_);
604 for(int i = 0; i < numberOfInputs_; i++) {
605 min_diag_price[i] = 0;
606 min_price[i] = 0;
607 }
608
609 int min_points_to_add = std::numeric_limits<int>::max();
610 this->enrichCurrentBidderDiagrams(2 * max_persistence, min_persistence,
611 min_diag_price, min_price,
612 min_points_to_add, false);
613
614 bool converged = false;
615 bool finished = false;
616 double total_cost;
617
618 while(!finished) {
619 Timer tm;
620
621 std::pair<std::unique_ptr<KDT>, std::vector<KDT *>> pair;
622 bool use_kdt = false;
623 // If the barycenter is empty, do not compute the kdt (or it will crash :/)
624 // TODO Fix KDTree to handle empty inputs...
625 if(!barycenter_goods_[0].empty()) {
626 pair = this->getKDTree();
627 use_kdt = true;
628 }
629
630 std::vector<std::vector<MatchingType>> all_matchings(numberOfInputs_);
631 std::vector<int> sizes(numberOfInputs_);
632 for(int i = 0; i < numberOfInputs_; i++) {
633 sizes[i] = current_bidder_diagrams_[i].size();
634 }
635
636 total_cost = 0;
637
638 barycenter.clear();
639 for(size_t j = 0; j < barycenter_goods_[0].size(); j++) {
640 Good &g = barycenter_goods_[0].at(j);
641 barycenter.emplace_back(PersistencePair{CriticalVertex{0, nt1_, g.x_, {}},
642 CriticalVertex{0, nt2_, g.y_, {}},
643 diagramType_, true});
644 }
645
646 bool actual_distance = (numberOfInputs_ == 2);
647 runMatchingAuction(&total_cost, sizes, *pair.first, pair.second,
648 &min_diag_price, &all_matchings, use_kdt,
649 actual_distance);
650
651 this->printMsg("Barycenter cost : " + std::to_string(total_cost),
653
654 if(converged) {
655 finished = true;
656 }
657
658 if(!finished) {
659 updateBarycenter(all_matchings);
660
661 if(min_cost > total_cost) {
662 min_cost = total_cost;
663 last_min_cost_obtained = 0;
664 } else {
665 last_min_cost_obtained += 1;
666 }
667
668 converged = converged || last_min_cost_obtained > 1;
669 }
670
671 previous_matchings = std::move(all_matchings);
672 // END OF TIMER
673
674 for(size_t i = 0; i < barycenter_goods_.size(); ++i) {
675 for(size_t j = 0; j < barycenter_goods_[i].size(); ++j) {
676 barycenter_goods_[i].at(j).setPrice(0);
677 }
678 }
679 for(size_t i = 0; i < current_bidder_diagrams_.size(); ++i) {
680 for(size_t j = 0; j < current_bidder_diagrams_[i].size(); ++j) {
681 current_bidder_diagrams_[i].at(j).setDiagonalPrice(0);
682 }
683 }
684 for(int i = 0; i < numberOfInputs_; i++) {
685 min_diag_price[i] = 0;
686 min_price[i] = 0;
687 }
688 }
689 barycenter.resize(0);
690 for(size_t j = 0; j < barycenter_goods_[0].size(); j++) {
691 Good &g = barycenter_goods_[0].at(j);
692 barycenter.emplace_back(PersistencePair{CriticalVertex{0, nt1_, g.x_, {}},
693 CriticalVertex{0, nt2_, g.y_, {}},
694 diagramType_, true});
695 }
696
697 cost_ = total_cost;
698 std::vector<std::vector<MatchingType>> corrected_matchings
699 = correctMatchings(previous_matchings);
700 return corrected_matchings;
701}
702
704 double total_real_cost = 0;
705 std::vector<MatchingType> fake_matchings;
706 for(int i = 0; i < numberOfInputs_; i++) {
708 wasserstein_, geometrical_factor_, lambda_, 0.01, true);
709 GoodDiagram current_barycenter = barycenter_goods_[0];
710 BidderDiagram current_bidder_diagram = bidder_diagrams_[i];
711 auction.BuildAuctionDiagrams(current_bidder_diagram, current_barycenter);
712 double cost = auction.run(fake_matchings);
713 total_real_cost += cost * cost;
714 }
715 return sqrt(total_real_cost);
716}
717
718bool ttk::PDBarycenter::isPrecisionObjectiveMet(double precision_objective,
719 int mode) {
720 if(mode == 0) { // ABSOLUTE PRECISION
721 for(int i_input = 0; i_input < numberOfInputs_; i_input++) {
722 if(precision_[i_input] > precision_objective) {
723 return false;
724 }
725 }
726 } else if(mode == 1) { // AVERAGE PRECISION
727 double average_precision
728 = std::accumulate(precision_.begin(), precision_.end(), 0.0)
729 / numberOfInputs_;
730 if(average_precision > precision_objective) {
731 return false;
732 }
733 }
734 return true;
735}
void setDiagonalPrice(const double price)
void setPositionInAuction(const int pos)
void setPrice(const double price)
TTK KD-Tree.
Definition: KDTree.h:21
std::vector< std::vector< MatchingType > > executeAuctionBarycenter(DiagramType &barycenter)
double getEpsilon(double rho)
void runMatching(double *total_cost, double epsilon, std::vector< int > &sizes, KDT &kdt, std::vector< KDT * > &correspondence_kdt_map, std::vector< double > *min_diag_price, std::vector< double > *min_price, std::vector< std::vector< MatchingType > > *all_matchings, bool use_kdt, bool actual_distance)
std::pair< typename KDT::KDTreeRoot, typename KDT::KDTreeMap > KDTreePair
Definition: PDBarycenter.h:58
double getRho(double epsilon)
double computeRealCost()
double updateBarycenter(std::vector< std::vector< MatchingType > > &matchings)
void setInitialBarycenter(double min_persistence)
bool isPrecisionObjectiveMet(double, int)
KDTreePair getKDTree() const
std::vector< std::vector< MatchingType > > correctMatchings(std::vector< std::vector< MatchingType > > &previous_matchings)
void runMatchingAuction(double *total_cost, std::vector< int > &sizes, KDT &kdt, std::vector< KDT * > &correspondence_kdt_map, std::vector< double > *min_diag_price, std::vector< std::vector< MatchingType > > *all_matchings, bool use_kdt, bool actual_distance)
double getMaxPersistence()
double getMinimalPrice(int i)
double getLowestPersistence()
std::vector< std::vector< MatchingType > > execute(DiagramType &barycenter)
bool hasBarycenterConverged(std::vector< std::vector< MatchingType > > &matchings, std::vector< std::vector< MatchingType > > &previous_matchings)
double enrichCurrentBidderDiagrams(double previous_min_persistence, double min_persistence, std::vector< double > &initial_diagonal_prices, std::vector< double > &initial_off_diagonal_prices, int min_points_to_add, bool add_points_to_barycenter=true)
void SetCriticalCoordinates(const float coords_x, const float coords_y, const float coords_z)
std::array< float, 3 > GetCriticalCoordinates() const
void runAuctionRound(int &n_biddings, const int kdt_index=0)
double run(std::vector< MatchingType > &matchings, const int kdt_index=0)
double initLowerBoundCost(const int kdt_index=0)
double getMatchingsAndDistance(std::vector< MatchingType > &matchings, bool get_diagonal_matches=false)
void BuildAuctionDiagrams(const BidderDiagram &BD, const GoodDiagram &GD)
void initLowerBoundCostWeight(double delta_lim)
double getElapsedTime()
Definition: Timer.h:15
T1 pow(const T1 val, const T2 n)
Definition: Geometry.h:411
std::vector< Bidder > BidderDiagram
std::vector< PersistencePair > DiagramType
Persistence Diagram type as a vector of Persistence pairs.
std::tuple< int, int, double > MatchingType
Matching between two Persistence Diagram pairs.
std::vector< Good > GoodDiagram
printMsg(debug::output::GREEN+"                           "+debug::output::ENDCOLOR+debug::output::GREEN+"▒"+debug::output::ENDCOLOR+debug::output::GREEN+"▒▒▒▒▒▒▒▒▒▒▒▒▒░"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)