TTK
Loading...
Searching...
No Matches
PersistenceDiagramAuction.cpp
Go to the documentation of this file.
2
4 const int kdt_index) {
5 double max_price = getMaximalPrice();
6 double epsilon = epsilon_;
7 if(epsilon_ < 1e-6 * max_price) {
8 // Risks of floating point limits reached...
9 epsilon = 1e-6 * max_price;
10 }
11 while(unassignedBidders_.size() > 0) {
12 n_biddings++;
13 int pos = unassignedBidders_.front();
14 Bidder &b = this->bidders_[pos];
16
17 GoodDiagram &all_goods = b.isDiagonal() ? diagonal_goods_ : goods_;
18 Good &twin_good = b.id_ >= 0 ? diagonal_goods_[b.id_] : goods_[-b.id_ - 1];
19 // double eps = epsilon_*(1+0.05*n_biddings/bidders_.size());
20 int idx_reassigned;
21 if(b.isDiagonal()) {
22 if(use_kdt_) {
23 idx_reassigned = b.runDiagonalKDTBidding(
24 &all_goods, twin_good, wasserstein_, epsilon, geometricalFactor_,
26 } else {
27 idx_reassigned
28 = b.runDiagonalBidding(&all_goods, twin_good, wasserstein_, epsilon,
30 }
31 } else {
32 if(use_kdt_) {
33 // We can use the kd-tree to speed up the search
34 idx_reassigned
35 = b.runKDTBidding(&all_goods, twin_good, wasserstein_, epsilon,
36 geometricalFactor_, &kdt_, kdt_index);
37 } else {
38 idx_reassigned = b.runBidding(
39 &all_goods, twin_good, wasserstein_, epsilon, geometricalFactor_);
40 }
41 }
42 if(idx_reassigned >= 0) {
43 Bidder &reassigned = bidders_[idx_reassigned];
44 reassigned.resetProperty();
45 unassignedBidders_.push(idx_reassigned);
46 }
47 }
48}
49
51 double max_price = 0;
52 for(size_t i = 0; i < goods_.size(); ++i) {
53 Good &g = goods_[i];
54 double price = g.getPrice();
55 if(price > max_price) {
56 max_price = price;
57 }
58 }
59
60 for(size_t i = 0; i < diagonal_goods_.size(); ++i) {
61 Good &g = diagonal_goods_[i];
62 double price = g.getPrice();
63 if(price > max_price) {
64 max_price = price;
65 }
66 }
67
68 return max_price;
69}
70
72 std::vector<MatchingType> &matchings, bool get_diagonal_matches) {
73 double wassersteinDistance = 0;
74 for(size_t i = 0; i < bidders_.size(); i++) {
75 Bidder &b = bidders_[i];
76 if(!b.isDiagonal()) {
77 int good_id = b.getProperty().id_;
78 double cost;
79
80 if(good_id > -1) {
81 // good is not diagonal
82 cost = b.cost(b.getProperty(), wasserstein_, geometricalFactor_);
83 MatchingType t = std::make_tuple(i, good_id, cost);
84 matchings.emplace_back(t);
85 } else {
86 if(!b.getProperty().isDiagonal()) {
87 this->printWrn("both the bidder and the good are diagonal points");
88 }
89 cost = 2.0 * Geometry::pow(std::abs((b.y_ - b.x_) / 2.0), wasserstein_);
90 if(get_diagonal_matches) {
91 MatchingType t = std::make_tuple(i, good_id, cost);
92 matchings.emplace_back(t);
93 }
94 }
95 wassersteinDistance += cost;
96 } else {
97 // b is diagonal
98 const Good &g = b.getProperty();
99 double cost
100 = 2.0 * Geometry::pow(std::abs((g.y_ - g.x_) / 2.0), wasserstein_);
101 if(get_diagonal_matches) {
102 MatchingType t = std::make_tuple(b.id_, g.id_, cost);
103 matchings.emplace_back(t);
104 }
105 wassersteinDistance += cost;
106 }
107 }
108 return wassersteinDistance;
109}
110
112 lowerBoundCost_ = 0;
113 for(unsigned int i = 0; i < bidders_.size(); ++i) {
114 if(bidders_[i].isDiagonal())
115 continue;
116
117 // Get closest good
118 double bestCost = std::numeric_limits<double>::max();
119 std::vector<KDT *> neighbours;
120 std::vector<double> costs;
121 if(use_kdt_) {
122 std::array<double, 5> coordinates;
123 bidders_[i].GetKDTCoordinates(geometricalFactor_, coordinates);
124 kdt_.getKClosest(1, coordinates, neighbours, costs, kdt_index);
125 int bestIndex = neighbours[0]->id_;
126 bestCost
127 = bidders_[i].cost(goods_[bestIndex], wasserstein_, geometricalFactor_);
128 } else {
129 for(unsigned int j = 0; j < goods_.size(); ++j) {
130 double cost
131 = bidders_[i].cost(goods_[j], wasserstein_, geometricalFactor_);
132 if(cost < bestCost)
133 bestCost = cost;
134 }
135 }
136
137 // Compare with diagonal good
138 Good g{bidders_[i].x_, bidders_[i].y_, true, -bidders_[i].id_ - 1};
139 g.projectOnDiagonal();
140 double cost = bidders_[i].cost(g, wasserstein_, geometricalFactor_);
141 if(cost < bestCost)
142 bestCost = cost;
143
144 // Update lower bound
145 lowerBoundCost_ += bestCost;
146 }
147 return lowerBoundCost_;
148}
149
150double ttk::PersistenceDiagramAuction::run(std::vector<MatchingType> &matchings,
151 const int kdt_index) {
152 initLowerBoundCostWeight(delta_lim_);
153 initLowerBoundCost(kdt_index);
154 initializeEpsilon();
155 int n_biddings = 0;
156 double delta = 5;
157 while(delta > delta_lim_) {
158 epsilon_ /= 5;
159 this->buildUnassignedBidders();
160 this->reinitializeGoods();
161 this->runAuctionRound(n_biddings, kdt_index);
162 delta = this->getRelativePrecision();
163 }
164 double wassersteinDistance = this->getMatchingsAndDistance(matchings, true);
165 return wassersteinDistance;
166}
167
170 const int wasserstein,
171 const double geometricalFactor) const {
172
173 if(is_diagonal_ && g.isDiagonal()) {
174 return 0;
175 } else if(is_diagonal_) {
176 return geometricalFactor
177 * (2 * Geometry::pow(std::abs(g.y_ / 2 - g.x_ / 2), wasserstein))
178 + (1 - geometricalFactor) * getPairGeometricalLength(wasserstein);
179 } else if(g.isDiagonal()) {
180 return geometricalFactor
181 * (2 * Geometry::pow(std::abs(y_ / 2 - x_ / 2), wasserstein))
182 + (1 - geometricalFactor) * g.getPairGeometricalLength(wasserstein);
183 } else {
184 return geometricalFactor
185 * (Geometry::pow(std::abs(x_ - g.x_), wasserstein)
186 + Geometry::pow(std::abs(y_ - g.y_), wasserstein))
187 + (1 - geometricalFactor)
188 * (Geometry::pow(
189 std::abs(coords_[0] - g.coords_[0]), wasserstein)
191 std::abs(coords_[1] - g.coords_[1]), wasserstein)
193 std::abs(coords_[2] - g.coords_[2]), wasserstein));
194 }
195}
196
198 Good &twinGood,
199 int wasserstein,
200 double epsilon,
201 double geometricalFactor) {
202 double best_val = std::numeric_limits<double>::lowest();
203 double second_val = std::numeric_limits<double>::lowest();
204 Good *best_good{};
205 for(size_t i = 0; i < goods->size(); i++) {
206 Good &g = (*goods)[i];
207 double val = -this->cost(g, wasserstein, geometricalFactor);
208 val -= g.getPrice();
209 if(val > best_val) {
210 second_val = best_val;
211 best_val = val;
212 best_good = &g;
213 } else if(val > second_val) {
214 second_val = val;
215 }
216 }
217 // And now check for the corresponding twin bidder
218 Good &g = twinGood;
219 double val = -this->cost(g, wasserstein, geometricalFactor);
220 val -= g.getPrice();
221 if(val > best_val) {
222 second_val = best_val;
223 best_val = val;
224 best_good = &g;
225 } else if(val > second_val) {
226 second_val = val;
227 }
228
229 if(second_val == std::numeric_limits<double>::lowest()) {
230 // There is only one acceptable good for the bidder
231 second_val = best_val;
232 }
233 if(best_good == nullptr) {
234 return -1;
235 }
236 double old_price = best_good->getPrice();
237 double new_price = old_price + best_val - second_val + epsilon;
238 if(new_price > std::numeric_limits<double>::max() / 2) {
239 new_price = old_price + epsilon;
240 }
241 // Assign bidder to best_good
242 this->setProperty(*best_good);
243 this->setPricePaid(new_price);
244
245 // Assign best_good to bidder and unassign the previous owner of best_good
246 // if need be
247 int idx_reassigned = best_good->getOwner();
248 best_good->assign(this->position_in_auction_, new_price);
249 return idx_reassigned;
250}
251
253 GoodDiagram *goods,
254 Good &twinGood,
255 int wasserstein,
256 double epsilon,
257 double geometricalFactor,
258 std::priority_queue<std::pair<int, double>,
259 std::vector<std::pair<int, double>>,
260 Compare> &diagonal_queue) {
261
262 // First, find the lowest and second lowest weights for diagonal goods
263 // Take this time to update the weights in the priority queue of the goods
264 // tested. It is not necessary to update weights for all diagonal bidders in
265 // the pririty queue since weights can only increase and we are interested
266 // only in the lowest ones
267 bool updated_top_pair
268 = false; // Boolean which equals true iff the top pair in the priority
269 // queue is given the good price
270 bool non_empty_goods = goods->size() > 0;
271 std::pair<int, double> best_pair;
272 if(non_empty_goods) {
273 while(!updated_top_pair) {
274 std::pair<int, double> top_pair = diagonal_queue.top();
275 diagonal_queue.pop();
276
277 double queue_weight = top_pair.second;
278 const auto &good = (*goods)[top_pair.first];
279 if(good.getPrice() > queue_weight) {
280 // If the weight in the priority queue is not the good one, update
281 std::get<1>(top_pair) = good.getPrice();
282 diagonal_queue.push(top_pair);
283 } else {
284 updated_top_pair = true;
285 best_pair = top_pair;
286 }
287 }
288 }
289
290 std::pair<int, double> second_pair;
291 if(non_empty_goods && diagonal_queue.size() > 0) {
292 bool updated_second_pair = false;
293 while(!updated_second_pair) {
294 second_pair = diagonal_queue.top();
295 double queue_weight = second_pair.second;
296 const auto &good = (*goods)[second_pair.first];
297 if(good.getPrice() != queue_weight) {
298 // If the weight in the priority queue is not the good one, update it
299 diagonal_queue.pop();
300 std::get<1>(second_pair) = good.getPrice();
301 diagonal_queue.push(second_pair);
302 } else {
303 updated_second_pair = true;
304 }
305 }
306 }
307
308 Good *best_good{};
309 double best_val = 0;
310 double second_val = 0;
311 if(non_empty_goods) {
312 best_val = -best_pair.second;
313 second_val = diagonal_queue.size() > 0 ? -second_pair.second : best_val;
314 best_good = &(*goods)[best_pair.first];
315 }
316
317 // And now check for the corresponding twin bidder
318 bool is_twin = false;
319 Good &g = twinGood;
320 double val = -this->cost(g, wasserstein, geometricalFactor);
321 val -= g.getPrice();
322 if(non_empty_goods) {
323 if(val > best_val) {
324 second_val = best_val;
325 best_val = val;
326 best_good = &g;
327 is_twin = true;
328 } else if(val > second_val || diagonal_queue.empty()) {
329 second_val = val;
330 }
331 } else {
332 best_val = val;
333 best_good = &g;
334 is_twin = true;
335 second_val = best_val;
336 }
337
338 if(second_val == std::numeric_limits<double>::lowest()) {
339 // There is only one acceptable good for the bidder
340 second_val = best_val;
341 }
342 double old_price = best_good->getPrice();
343 double new_price = old_price + best_val - second_val + epsilon;
344 if(new_price > std::numeric_limits<double>::max() / 2) {
345 new_price = old_price + epsilon;
346 }
347
348 // Assign bidder to best_good
349 this->setProperty(*best_good);
350 this->setPricePaid(new_price);
351
352 // Assign best_good to bidder and unassign the previous owner of best_good
353 // if need be
354 int idx_reassigned = best_good->getOwner();
355 best_good->assign(this->position_in_auction_, new_price);
356 if(!is_twin) {
357 std::get<1>(best_pair) = new_price;
358 }
359 if(non_empty_goods) {
360 diagonal_queue.push(best_pair);
361 }
362 return idx_reassigned;
363}
364
366 GoodDiagram *goods,
367 Good &twinGood,
368 int wasserstein,
369 double epsilon,
370 double geometricalFactor,
371 std::vector<KDT *> &correspondence_kdt_map,
372 std::priority_queue<std::pair<int, double>,
373 std::vector<std::pair<int, double>>,
374 Compare> &diagonal_queue,
375 const int kdt_index) {
377
378 // First, find the lowest and second lowest weights for diagonal goods
379 // Take this time to update the weights in the priority queue of the goods
380 // tested. It is not necessary to update weights for all diagonal bidders in
381 // the pririty queue since weights can only increase and we are interested
382 // only in the lowest ones
383 bool updated_top_pair
384 = false; // Boolean which equals true iff the top pair in the priority
385 // queue is given the good price
386 bool non_empty_goods = goods->size() > 0;
387 std::pair<int, double> best_pair;
388 if(non_empty_goods) {
389 while(!updated_top_pair) {
390 std::pair<int, double> top_pair = diagonal_queue.top();
391 double queue_weight = top_pair.second;
392 const auto &good = (*goods)[top_pair.first];
393
394 diagonal_queue.pop();
395 if(good.getPrice() > queue_weight) {
396 // If the weight in the priority queue is not the good one, update it
397 std::get<1>(top_pair) = good.getPrice();
398 diagonal_queue.push(top_pair);
399 } else {
400 updated_top_pair = true;
401 best_pair = top_pair;
402 }
403 }
404 }
405 std::pair<int, double> second_pair;
406 if(non_empty_goods && diagonal_queue.size() > 0) {
407 bool updated_second_pair = false;
408 while(!updated_second_pair) {
409 second_pair = diagonal_queue.top();
410 double queue_weight = second_pair.second;
411 const auto &good = (*goods)[second_pair.first];
412 if(good.getPrice() > queue_weight) {
413 // If the weight in the priority queue is not the good one, update it
414 diagonal_queue.pop();
415 std::get<1>(second_pair) = good.getPrice();
416 diagonal_queue.push(second_pair);
417 } else {
418 updated_second_pair = true;
419 }
420 }
421 }
422
423 Good *best_good{};
424 double best_val = 0;
425 double second_val = 0;
426 if(non_empty_goods) {
427 best_val = -best_pair.second;
428 second_val = diagonal_queue.size() > 0 ? -second_pair.second : best_val;
429 best_good = &(*goods)[best_pair.first];
430 }
431
432 // And now check for the corresponding twin bidder
433 bool is_twin = false;
434 Good &g = twinGood;
435 double val = -this->cost(g, wasserstein, geometricalFactor);
436 val -= g.getPrice();
437 if(non_empty_goods) {
438 if(val > best_val) {
439 second_val = best_val;
440 best_val = val;
441 best_good = &g;
442 is_twin = true;
443 } else if(val > second_val || diagonal_queue.empty()) {
444 second_val = val;
445 }
446 } else {
447 best_val = val;
448 best_good = &g;
449 is_twin = true;
450 second_val = best_val;
451 }
452
453 if(second_val == std::numeric_limits<double>::lowest()) {
454 // There is only one acceptable good for the bidder
455 second_val = best_val;
456 }
457 double old_price = best_good->getPrice();
458 double new_price = old_price + best_val - second_val + epsilon;
459 if(new_price > std::numeric_limits<double>::max() / 2) {
460 new_price = old_price + epsilon;
461 }
462
463 // Assign bidder to best_good
464 this->setProperty(*best_good);
465 this->setPricePaid(new_price);
466
467 // Assign best_good to bidder and unassign the previous owner of best_good
468 // if need be
469 int idx_reassigned = best_good->getOwner();
470 best_good->assign(this->position_in_auction_, new_price);
471 if(is_twin) {
472 // Update weight in KDTree if the closest good is in it
473 correspondence_kdt_map[best_good->id_]->updateWeight(new_price, kdt_index);
474 if(non_empty_goods) {
475 diagonal_queue.push(best_pair);
476 }
477 } else {
478 if(non_empty_goods) {
479 std::get<1>(best_pair) = new_price;
480 diagonal_queue.push(best_pair);
481 }
482 }
483 return idx_reassigned;
484}
485
487 Good &twinGood,
488 int wasserstein,
489 double epsilon,
490 double geometricalFactor,
491 KDT *kdt,
492 const int kdt_index) {
493
495 std::vector<KDT *> neighbours;
496 std::vector<double> costs;
497
498 std::array<double, 5> coordinates;
499 GetKDTCoordinates(geometricalFactor, coordinates);
500
501 kdt->getKClosest(2, coordinates, neighbours, costs, kdt_index);
502 double best_val, second_val;
503 KDT *closest_kdt;
504 Good *best_good{};
505 if(costs.size() == 2) {
506 std::array<int, 2> idx{0, 1};
507 std::sort(idx.begin(), idx.end(),
508 [&costs](int &a, int &b) { return costs[a] < costs[b]; });
509
510 closest_kdt = neighbours[idx[0]];
511 best_good = &(*goods)[closest_kdt->id_];
512 // Value is defined as the opposite of cost (each bidder aims at
513 // maximizing it)
514 best_val = -costs[idx[0]];
515 second_val = -costs[idx[1]];
516 } else {
517 // If the kdtree contains only one point
518 closest_kdt = neighbours[0];
519 best_good = &(*goods)[closest_kdt->id_];
520 best_val = -costs[0];
521 second_val = best_val;
522 }
523 // And now check for the corresponding twin bidder
524 bool twin_chosen = false;
525 Good &g = twinGood;
526 double val = -this->cost(g, wasserstein, geometricalFactor);
527 val -= g.getPrice();
528 if(val > best_val) {
529 second_val = best_val;
530 best_val = val;
531 best_good = &g;
532 twin_chosen = true;
533 } else if(val > second_val) {
534 second_val = val;
535 }
536
537 if(second_val == std::numeric_limits<double>::lowest()) {
538 // There is only one acceptable good for the bidder
539 second_val = best_val;
540 }
541 double old_price = best_good->getPrice();
542 double new_price = old_price + best_val - second_val + epsilon;
543 if(new_price > std::numeric_limits<double>::max() / 2) {
544 new_price = old_price + epsilon;
545 }
546 // Assign bidder to best_good
547 this->setProperty(*best_good);
548 this->setPricePaid(new_price);
549 // Assign best_good to bidder and unassign the previous owner of best_good
550 // if need be
551 int idx_reassigned = best_good->getOwner();
552
553 best_good->assign(this->position_in_auction_, new_price);
554 // Update the price in the KDTree
555 if(!twin_chosen) {
556 closest_kdt->updateWeight(new_price, kdt_index);
557 }
558 return idx_reassigned;
559}
int runKDTBidding(GoodDiagram *goods, Good &diagonalGood, int wasserstein, double epsilon, double geometricalFactor, KDT *kdt, const int kdt_index=0)
int runBidding(GoodDiagram *goods, Good &diagonalGood, int wasserstein, double epsilon, double geometricalFactor)
int runDiagonalKDTBidding(GoodDiagram *goods, Good &diagonalGood, int wasserstein, double epsilon, double geometricalFactor, std::vector< KDT * > &correspondence_kdt_map, std::priority_queue< std::pair< int, double >, std::vector< std::pair< int, double > >, Compare > &diagonal_queue, const int kdt_index=0)
int runDiagonalBidding(GoodDiagram *goods, Good &twinGood, int wasserstein, double epsilon, double geometricalFactor, std::priority_queue< std::pair< int, double >, std::vector< std::pair< int, double > >, Compare > &diagonal_queue)
const Good & getProperty() const
TTK KD-Tree.
Definition: KDTree.h:21
void updateWeight(const dataType new_weight, const int weight_index=0)
Definition: KDTree.h:78
void getKClosest(const unsigned int k, const Container &coordinates, KDTreeMap &neighbours, std::vector< dataType > &costs, const int weight_index=0)
Definition: KDTree.h:335
int id_
Definition: KDTree.h:44
double cost(const PersistenceDiagramAuctionActor &g, const int wasserstein, const double geometricalFactor) const
double getPairGeometricalLength(const int wasserstein) const
void runAuctionRound(int &n_biddings, const int kdt_index=0)
std::vector< KDT * > & correspondence_kdt_map_
std::priority_queue< std::pair< int, double >, std::vector< std::pair< int, double > >, Compare > diagonal_queue_
double initLowerBoundCost(const int kdt_index=0)
double getMatchingsAndDistance(std::vector< MatchingType > &matchings, bool get_diagonal_matches=false)
T1 pow(const T1 val, const T2 n)
Definition: Geometry.h:411
std::tuple< int, int, double > MatchingType
Matching between two Persistence Diagram pairs.
std::vector< Good > GoodDiagram