TTK
Loading...
Searching...
No Matches
PDClustering.cpp
Go to the documentation of this file.
1
14
15#include <PDClustering.h>
16
17#include <algorithm>
18#include <cmath>
19#include <cstdlib> /* srand, rand */
20#include <iostream>
21#include <iterator>
22#include <random>
23
25 std::vector<DiagramType> &final_centroids,
26 std::vector<std::vector<std::vector<std::vector<MatchingType>>>>
27 &all_matchings_per_type_and_cluster) {
28
29 all_matchings_per_type_and_cluster.resize(k_);
30 for(int c = 0; c < k_; c++) {
31 all_matchings_per_type_and_cluster[c].resize(3);
32 for(int i = 0; i < 3; i++) {
33 all_matchings_per_type_and_cluster[c][i].resize(numberOfInputs_);
34 }
35 }
36 int matchings_only = false;
37 Timer tm;
38 {
39 // PARTICULARITIES FOR THE CASE OF ONE UNIQUE CLUSTER
40 if(k_ <= 1) {
41 use_accelerated_ = false;
42 use_kmeanspp_ = false;
44 use_progressive_ = false;
45 deterministic_ = true;
46 matchings_only = true;
47 time_limit_ = 99999999999;
48 }
49 }
50
51 std::vector<bool *> current_prec;
52 current_prec.emplace_back(&precision_min_);
53 current_prec.emplace_back(&precision_sad_);
54 current_prec.emplace_back(&precision_max_);
55
56 std::vector<bool *> current_dos;
57 current_dos.emplace_back(&do_min_);
58 current_dos.emplace_back(&do_sad_);
59 current_dos.emplace_back(&do_max_);
60 bool converged = false;
61 std::vector<bool> diagrams_complete(3);
62 for(int c = 0; c < 3; c++) {
63 diagrams_complete[c] = (!use_progressive_) || (!original_dos[c]);
64 }
65 n_iterations_ = 0;
66 double total_time = 0;
67
68 // double cost = std::numeric_limits<double>::max();
70 cost_ = std::numeric_limits<double>::max();
71 double min_cost_min = std::numeric_limits<double>::max();
72 double min_cost_max = std::numeric_limits<double>::max();
73 double min_cost_sad = std::numeric_limits<double>::max();
74 // double last_min_cost_obtained = -1;
75 double last_min_cost_obtained_min = -1;
76 double last_min_cost_obtained_sad = -1;
77 double last_min_cost_obtained_max = -1;
78 std::vector<double> epsilon0(3);
79 std::vector<double> epsilon_candidate(3);
80 std::vector<double> rho(3);
81
82 // Getting current diagrams (with only at most min_points_to_add points)
83 std::vector<double> max_persistence(3);
84 std::vector<double> lowest_persistence(3);
85 std::vector<double> min_persistence(3);
86
87 const auto square = [](const double a) { return a * a; };
88
89 for(int i_crit = 0; i_crit < 3; i_crit++) {
90 max_persistence[i_crit] = 2 * getMostPersistent(i_crit);
91 lowest_persistence[i_crit] = getLessPersistent(i_crit);
92 min_persistence[i_crit] = 0;
93 epsilon_[i_crit] = square(0.5 * max_persistence[i_crit])
94 / 8.; // max_persistence actually holds 2 times the
95 // highest persistence
96 epsilon0[i_crit] = epsilon_[i_crit];
97 }
98
99 std::vector<int> min_points_to_add(3);
100 min_points_to_add[0] = 10;
101 min_points_to_add[1] = 10;
102 min_points_to_add[2] = 10;
103
104 if(use_progressive_) {
105 } else {
106 min_points_to_add[0] = std::numeric_limits<int>::max();
107 min_points_to_add[1] = std::numeric_limits<int>::max();
108 min_points_to_add[2] = std::numeric_limits<int>::max();
109 }
110 std::vector<std::vector<double>> min_diag_price(3);
111 std::vector<std::vector<double>> min_off_diag_price(3);
112 for(int c = 0; c < 3; ++c) {
113 for(int i = 0; i < numberOfInputs_; i++) {
114 min_diag_price[c].emplace_back(0);
115 min_off_diag_price[c].emplace_back(0);
116 }
117 }
118 min_persistence = enrichCurrentBidderDiagrams(
119 max_persistence, min_persistence, min_diag_price, min_off_diag_price,
120 min_points_to_add, false, true);
121
122 for(int c = 0; c < 3; c++) {
123 if(min_persistence[c] <= lowest_persistence[c]) {
124 diagrams_complete[c] = true;
125 } // max_persistence actually holds 2 times the highest persistence
126 }
127 bool all_diagrams_complete
128 = diagrams_complete[0] && diagrams_complete[1] && diagrams_complete[2];
129 if(all_diagrams_complete) {
130 use_progressive_ = false;
131 }
132
133 // Initializing centroids and clusters
134 if(use_kmeanspp_) {
136 } else {
138 }
140 if(use_accelerated_) {
144 } else {
147 }
148 if(debugLevel_ > 3 && k_ > 1) {
149 printMsg("Initial Clustering: ");
151 }
152 initializeBarycenterComputers(min_persistence);
153 while(!converged || (!all_diagrams_complete && use_progressive_)) {
154 Timer t_inside;
155 {
157
158 for(int i_crit = 0; i_crit < 3; i_crit++) {
159 if(*(current_dos[i_crit])) {
160 rho[i_crit] = min_persistence[i_crit] > 0
161 ? std::sqrt(8.0 * epsilon_[i_crit])
162 : -1;
163 }
164 }
165
167
168 do_min_ = do_min_ && (min_persistence[0] > rho[0]);
169 do_sad_ = do_sad_ && (min_persistence[1] > rho[1]);
170 do_max_ = do_max_ && (min_persistence[2] > rho[2]);
171
172 for(int i_crit = 0; i_crit < 3; i_crit++) {
173 if(*(current_dos[i_crit])) {
174 epsilon_candidate[i_crit] = square(min_persistence[i_crit]) / 8.;
175 if(epsilon_candidate[i_crit] > epsilon_[i_crit]) {
176 // Should always be the case except if min_persistence is
177 // equal to zero
178 epsilon_[i_crit] = epsilon_candidate[i_crit];
179 }
180 }
181 }
182
183 if(epsilon_[0] < 5e-5) {
184 // Add all remaining points for final convergence.
185 // rho[0] = 0;
186 min_persistence[0] = 0;
187 min_points_to_add[0] = std::numeric_limits<int>::max();
188 }
189 if(epsilon_[1] < 5e-5) {
190 // Add all remaining points for final convergence.
191 // rho[1] = 0;
192 min_persistence[1] = 0;
193 min_points_to_add[1] = std::numeric_limits<int>::max();
194 }
195 if(epsilon_[2] < 5e-5) {
196 // Add all remaining points for final convergence.
197 // rho[2] = 0;
198 min_persistence[2] = 0;
199 min_points_to_add[2] = std::numeric_limits<int>::max();
200 }
201
202 if(do_min_ || do_sad_ || do_max_) {
203 min_persistence = enrichCurrentBidderDiagrams(
204 min_persistence, rho, min_diag_price, min_off_diag_price,
205 min_points_to_add, true, false);
206 }
208
209 for(int i_crit = 0; i_crit < 3; i_crit++) {
210 if(*(current_dos[i_crit])) {
211 if(min_persistence[i_crit] <= lowest_persistence[i_crit]) {
212 diagrams_complete[i_crit] = true;
213 }
214 }
215 }
216
217 if(diagrams_complete[0] && diagrams_complete[1]
218 && diagrams_complete[2]) {
219 use_progressive_ = false;
220 all_diagrams_complete = true;
221 }
222
224 }
225 std::vector<double> max_shift_vec = updateCentroidsPosition(
226 &min_off_diag_price, &min_diag_price,
227 all_matchings_per_type_and_cluster, matchings_only);
228 if(do_min_ && !UseDeltaLim_) {
229 precision_min_ = (epsilon_[0] < epsilon0[0] / 500.);
230 }
231 if(do_sad_ && !UseDeltaLim_) {
232 precision_sad_ = (epsilon_[1] < epsilon0[1] / 500.);
233 }
234 if(do_max_ && !UseDeltaLim_) {
235 precision_max_ = (epsilon_[2] < epsilon0[2] / 500.);
236 }
237
238 for(int i_crit = 0; i_crit < 3; i_crit++) {
239 if(*(current_dos[i_crit]) /*&& (!*(current_prec[i_crit]) || !diagrams_complete[i_crit] ) */) {
240 epsilon_candidate[i_crit] = std::min(
241 std::max(max_shift_vec[i_crit] / 8., epsilon_[i_crit] / 5.),
242 epsilon0[i_crit] / square(n_iterations_));
243
244 if((epsilon_candidate[i_crit] < epsilon_[i_crit]
245 && !diagrams_complete[i_crit])
246 || diagrams_complete[i_crit]) {
247 epsilon_[i_crit] = epsilon_candidate[i_crit];
248 } else {
249 epsilon_[i_crit] *= 0.95;
250 }
251 }
252 }
253
254 if(epsilon_[0] < epsilon_min_ /*&& diagrams_complete[0]*/) {
255 this->printMsg("[min barycenter] epsilon under minimal value ",
257 do_min_ = false;
259 diagrams_complete[0] = true;
260 }
261 if(epsilon_[1] < epsilon_min_ /*&& diagrams_complete[1]*/) {
262 this->printMsg("[sad barycenter] epsilon under minimal value ",
264 do_sad_ = false;
266 diagrams_complete[1] = true;
267 }
268 if(epsilon_[2] < epsilon_min_ /*&& diagrams_complete[2]*/) {
269 this->printWrn("[max barycenter] epsilon under minimal value ");
270 do_max_ = false;
272 diagrams_complete[2] = true;
273 }
274
275 if(diagrams_complete[0] && diagrams_complete[1]
276 && diagrams_complete[2]) {
277 use_progressive_ = false;
278 all_diagrams_complete = true;
279 }
280 if(use_accelerated_) {
282 } else {
283 // updateClusters();
284 }
285
288 bool precision_criterion_reached = precision_criterion_;
289
290 this->printMsg("Iteration " + std::to_string(n_iterations_)
291 + " epsilon " + std::to_string(epsilon_[0]) + " "
292 + std::to_string(epsilon_[1]) + " "
293 + std::to_string(epsilon_[2]),
295 this->printMsg(" complete " + std::to_string(diagrams_complete[0]) + " "
296 + std::to_string(diagrams_complete[1]) + " "
297 + std::to_string(diagrams_complete[2]),
299 this->printMsg(" precision " + std::to_string(precision_min_) + " "
300 + std::to_string(precision_sad_) + " "
301 + std::to_string(precision_max_),
303 this->printMsg(" cost " + std::to_string(cost_min_) + " "
304 + std::to_string(cost_sad_) + " "
305 + std::to_string(cost_max_),
307
308 if(cost_min_ < min_cost_min && n_iterations_ > 2
309 && diagrams_complete[0] /*&& precision_min_*/) {
310 min_cost_min = cost_min_;
311 last_min_cost_obtained_min = 0;
312 } else if(n_iterations_ > 2 && precision_min_ && diagrams_complete[0]) {
313 last_min_cost_obtained_min += 1;
314 if(last_min_cost_obtained_min > 1) {
315 do_min_ = false;
316 }
317 }
318
319 if(cost_sad_ < min_cost_sad && n_iterations_ > 2
320 && diagrams_complete[1] /*&& precision_sad_*/) {
321 min_cost_sad = cost_sad_;
322 last_min_cost_obtained_sad = 0;
323 } else if(n_iterations_ > 2 && precision_sad_ && diagrams_complete[1]) {
324 last_min_cost_obtained_sad += 1;
325 if(last_min_cost_obtained_sad > 1 && diagrams_complete[1]) {
326 do_sad_ = false;
327 }
328 }
329
330 if(cost_max_ < min_cost_max && n_iterations_ > 2
331 && diagrams_complete[2] /*&& precision_max_*/) {
332 min_cost_max = cost_max_;
333 last_min_cost_obtained_max = 0;
334 } else if(n_iterations_ > 2 && precision_max_ && diagrams_complete[2]) {
335 last_min_cost_obtained_max += 1;
336 if(last_min_cost_obtained_max > 1 && diagrams_complete[2]) {
337 do_max_ = false;
338 }
339 }
340
341 if(debugLevel_ > 5) {
342 this->printMsg("Clustering result:", debug::Priority::DETAIL);
344 }
345 converged = converged
346 || (all_diagrams_complete && !do_min_ && !do_sad_
347 && !do_max_ && (precision_criterion_reached))
348 || (precision_criterion_reached && UseDeltaLim_
349 && numberOfInputs_ == 2);
350 }
351
352 total_time
353 += t_inside.getElapsedTime(); // - t_real_cost.getElapsedTime();
354
355 if(total_time + t_inside.getElapsedTime() > 0.9 * time_limit_) {
356 min_cost_min = cost_min_;
357 min_cost_sad = cost_sad_;
358 min_cost_max = cost_max_;
359 converged = true;
360 }
361 if(total_time > 0.1 * time_limit_) {
362 all_diagrams_complete = true;
363 diagrams_complete[0] = true;
364 diagrams_complete[1] = true;
365 diagrams_complete[2] = true;
366 use_progressive_ = false;
367 }
368 if(debugLevel_ > 4) {
369 std::cout << "== Iteration " << n_iterations_
370 << " == complete : " << all_diagrams_complete
371 << " , progressive : " << use_progressive_
372 << " , converged : " << converged << std::endl;
373 }
374 }
376
377 // display results
378 std::vector<std::vector<std::string>> rows{
379 {" Min-saddle cost", std::to_string(cost_min_)},
380 {" Saddle-saddle cost", std::to_string(cost_sad_)},
381 {" Saddle-max cost", std::to_string(cost_max_)},
382 {matchings_only ? "Wasserstein Distance" : "Final Cost",
383 std::to_string(cost_min_ + cost_sad_ + cost_max_)},
384 };
385 this->printMsg(rows);
386
387 if(!use_progressive_ && k_ > 1) {
388 clustering_ = old_clustering_; // reverting to last clustering
389 }
390 invertClusters(); // this is to pass the old inverse clustering to the VTK
391 // wrapper
392 if(k_ > 1) {
393 this->printMsg("Clustering result:");
395 }
396 } // End of timer
397
398 // CORRECT MATCHINGS :
399 // correctMatchings(all_matchings);
400 if(matchings_only) {
401 computeBarycenterForTwoGlobal(all_matchings_per_type_and_cluster);
402 }
403 correctMatchings(all_matchings_per_type_and_cluster);
404 // Filling the final centroids for output
405
406 final_centroids.resize(k_);
407 centroids_sizes_.resize(k_);
408 for(int c = 0; c < k_; c++) {
409 centroids_sizes_[c].resize(3);
410 if(do_min_)
411 centroids_sizes_[c][0] = centroids_min_[c].size();
412 if(do_sad_)
413 centroids_sizes_[c][1] = centroids_saddle_[c].size();
414 if(do_max_)
415 centroids_sizes_[c][2] = centroids_max_[c].size();
416 }
417
418 for(int c = 0; c < k_; ++c) {
419 // if NumberOfClusters > 1, the global pair was duplicated
420 // and needs to be removed from the min-saddle problem
421 // It is the first pair.
422 int removeFirstPairMin
423 = (k_ > 1 and original_dos[0] and original_dos[2]) ? 1 : 0;
424 int addedFirstPairMax = 0;
425 int addedFirstPairMin = removeFirstPairMin;
426
427 // min-max Pair
428 if(removeFirstPairMin or (!do_min_ and do_max_)) {
429 Good &g = centroids_max_[c].at(0);
430 const auto &critCoords = g.GetCriticalCoordinates();
431 final_centroids[c].emplace_back(PersistencePair{
433 CriticalVertex{0, CriticalType::Local_maximum, g.y_, critCoords}, 0,
434 false});
435 addedFirstPairMax = 1;
436 } else if(do_min_) {
437 Good &g = centroids_min_[c].at(0);
438 const auto &critCoords = g.GetCriticalCoordinates();
439 final_centroids[c].emplace_back(PersistencePair{
441 CriticalVertex{0, CriticalType::Local_maximum, g.y_, critCoords}, 0,
442 false});
443 addedFirstPairMin = 1;
444 }
445
446 if(do_min_) {
447 for(size_t i = addedFirstPairMin; i < centroids_min_[c].size(); ++i) {
448 Good &g = centroids_min_[c].at(i);
449 const auto &critCoords = g.GetCriticalCoordinates();
450 final_centroids[c].emplace_back(PersistencePair{
452 CriticalVertex{0, CriticalType::Saddle1, g.y_, critCoords}, 0, true});
453 if(g.getPersistence() > 1000) {
454 this->printMsg("Found a abnormally high persistence in min diagram",
456 }
457 }
458 }
459
460 if(do_sad_) {
461 for(size_t i = 0; i < centroids_saddle_[c].size(); ++i) {
462 Good &g = centroids_saddle_[c].at(i);
463 const auto &critCoords = g.GetCriticalCoordinates();
464 final_centroids[c].emplace_back(PersistencePair{
465 CriticalVertex{0, CriticalType::Saddle1, g.x_, critCoords},
466 CriticalVertex{0, CriticalType::Saddle2, g.y_, critCoords}, 1, true});
467 if(g.getPersistence() > 1000) {
468 this->printMsg("Found a abnormally high persistence in sad diagram",
470 }
471 }
472 }
473
474 if(do_max_) {
475 for(size_t i = addedFirstPairMax; i < centroids_max_[c].size(); ++i) {
476 Good &g = centroids_max_[c].at(i);
477 const auto &critCoords = g.GetCriticalCoordinates();
478 ttk::CriticalType saddle_type;
479 if(do_sad_)
480 saddle_type = ttk::CriticalType::Saddle2;
481 else
482 saddle_type = ttk::CriticalType::Saddle1;
483
484 final_centroids[c].emplace_back(PersistencePair{
485 CriticalVertex{0, saddle_type, g.x_, critCoords},
486 CriticalVertex{0, CriticalType::Local_maximum, g.y_, critCoords}, 2,
487 true});
488
489 if(g.getPersistence() > 1000) {
490 this->printMsg("Found a abnormally high persistence in min diagram",
492 }
493 }
494 }
495 }
496
497 if(distanceWritingOptions_ == 1) {
499 } else if(distanceWritingOptions_ == 2) {
501 }
502
503 return inv_clustering_;
504}
505
507 std::vector<std::vector<std::vector<std::vector<MatchingType>>>>
508 &previous_matchings) {
509 for(int c = 0; c < k_; c++) {
510 for(size_t i = 0; i < clustering_[c].size(); i++) {
511 int diagram_id = clustering_[c][i];
512 if(original_dos[0]) {
513 // 1. Invert the current_bidder_ids_ vector
514 std::vector<int> new_to_old_id(
515 current_bidder_diagrams_min_[diagram_id].size(), -1);
516 for(size_t j = 0; j < current_bidder_ids_min_[diagram_id].size(); j++) {
517 int new_id = current_bidder_ids_min_[diagram_id][j];
518 if(new_id >= 0) {
519 new_to_old_id[new_id] = j;
520 }
521 }
522 // 2. Reconstruct the matchings
523 std::vector<MatchingType> matchings_diagram_i;
524 for(size_t j = 0; j < previous_matchings[c][0][i].size(); j++) {
525 MatchingType m = previous_matchings[c][0][i][j];
526 int new_id = std::get<0>(m);
527 if(new_id >= 0 && std::get<1>(m) >= 0) {
528 std::get<0>(m) = new_to_old_id[new_id];
529 matchings_diagram_i.emplace_back(m);
530 } else if(std::get<1>(m)
531 >= 0) { // new_id < 0 corresponds to a diagonal matching
532 std::get<0>(m) = -1;
533 matchings_diagram_i.emplace_back(m);
534 }
535 }
536 previous_matchings[c][0][i].resize(matchings_diagram_i.size());
537 previous_matchings[c][0][i] = matchings_diagram_i;
538 }
539 if(original_dos[1]) {
540 // 1. Invert the current_bidder_ids_ vector
541 std::vector<int> new_to_old_id(
542 current_bidder_diagrams_saddle_[diagram_id].size());
543 for(size_t j = 0; j < current_bidder_ids_sad_[diagram_id].size(); j++) {
544 int new_id = current_bidder_ids_sad_[diagram_id][j];
545 if(new_id >= 0) {
546 new_to_old_id[new_id] = j;
547 }
548 }
549 // 2. Reconstruct the matchings
550 int zero_done = 0;
551 std::vector<MatchingType> matchings_diagram_i;
552 for(size_t j = 0; j < previous_matchings[c][1][i].size(); j++) {
553 MatchingType m = previous_matchings[c][1][i][j];
554 int new_id = std::get<0>(m);
555 if(new_id >= 0 && std::get<1>(m) >= 0) {
556 int old_id = new_to_old_id[new_id];
557 if(old_id > 0) {
558 std::get<0>(m) = old_id;
559 matchings_diagram_i.emplace_back(m);
560 } else {
561 if(!zero_done) {
562 zero_done = 1;
563 std::get<0>(m) = old_id;
564 matchings_diagram_i.emplace_back(m);
565 }
566 }
567 } else if(std::get<1>(m)
568 >= 0) { // new_id < 0 corresponds to a diagonal matching
569 std::get<0>(m) = -1;
570 matchings_diagram_i.emplace_back(m);
571 }
572 }
573 previous_matchings[c][1][i].resize(matchings_diagram_i.size());
574 previous_matchings[c][1][i] = matchings_diagram_i;
575 }
576 if(original_dos[2]) {
577 // 1. Invert the current_bidder_ids_ vector
578 std::vector<int> new_to_old_id(
579 current_bidder_diagrams_max_[diagram_id].size());
580 for(size_t j = 0; j < current_bidder_ids_max_[diagram_id].size(); j++) {
581 int new_id = current_bidder_ids_max_[diagram_id][j];
582 if(new_id >= 0) {
583 new_to_old_id[new_id] = j;
584 }
585 }
586 // 2. Reconstruct the matchings
587 int zero_done = 0;
588 std::vector<MatchingType> matchings_diagram_i;
589 for(size_t j = 0; j < previous_matchings[c][2][i].size(); j++) {
590 MatchingType m = previous_matchings[c][2][i][j];
591 int new_id = std::get<0>(m);
592 if(new_id >= 0 && std::get<1>(m) >= 0) {
593 int old_id = new_to_old_id[new_id];
594 if(old_id > 0) {
595 std::get<0>(m) = old_id;
596 matchings_diagram_i.emplace_back(m);
597 } else {
598 if(!zero_done) {
599 zero_done = 1;
600 std::get<0>(m) = old_id;
601 matchings_diagram_i.emplace_back(m);
602 }
603 }
604 } else if(std::get<1>(m)
605 >= 0) { // new_id < 0 corresponds to a diagonal matching
606 std::get<0>(m) = -1;
607 matchings_diagram_i.emplace_back(m);
608 }
609 }
610 previous_matchings[c][2][i].resize(matchings_diagram_i.size());
611 previous_matchings[c][2][i] = matchings_diagram_i;
612 }
613 }
614 }
615}
616
618 std::vector<std::vector<std::vector<MatchingType>>> matchings) {
619 std::cout << "\n MATCHINGS : " << std::endl;
620 for(int d = 0; d < 3; d++) {
621 if(original_dos[d]) {
622 std::cout << "\n Diagram type : " << d << std::endl;
623 for(size_t i = 0; i < matchings[d].size(); i++) {
624 std::cout << " diagram " << i << " : ";
625 for(size_t j = 0; j < matchings[d][i].size(); j++) {
626 std::cout << std::get<0>(matchings[d][i][j]) << " ";
627 std::cout << std::get<1>(matchings[d][i][j]) << " ";
628 std::cout << std::get<2>(matchings[d][i][j]) << " | ";
629 }
630 std::cout << "\n";
631 }
632 }
633 }
634}
635
636std::vector<std::vector<int>> ttk::PDClustering::get_centroids_sizes() {
637 return centroids_sizes_;
638}
639
641 double max_persistence = 0;
642
643 if(do_min_ && (type == -1 || type == 0)) {
644 for(size_t i = 0; i < bidder_diagrams_min_.size(); ++i) {
645 for(size_t j = 0; j < bidder_diagrams_min_[i].size(); ++j) {
646 Bidder b = bidder_diagrams_min_[i].at(j);
647 double persistence = b.getPersistence();
648 if(persistence > max_persistence) {
649 max_persistence = persistence;
650 }
651 }
652 }
653 }
654
655 if(do_sad_ && (type == -1 || type == 1)) {
656 for(size_t i = 0; i < bidder_diagrams_saddle_.size(); ++i) {
657 for(size_t j = 0; j < bidder_diagrams_saddle_[i].size(); ++j) {
658 Bidder b = bidder_diagrams_saddle_[i].at(j);
659 double persistence = b.getPersistence();
660 if(persistence > max_persistence) {
661 max_persistence = persistence;
662 }
663 }
664 }
665 }
666
667 if(do_max_ && (type == -1 || type == 2)) {
668 for(size_t i = 0; i < bidder_diagrams_max_.size(); ++i) {
669 for(size_t j = 0; j < bidder_diagrams_max_[i].size(); ++j) {
670 Bidder b = bidder_diagrams_max_[i].at(j);
671 double persistence = b.getPersistence();
672 if(persistence > max_persistence) {
673 max_persistence = persistence;
674 }
675 }
676 }
677 }
678 return max_persistence;
679}
680
682 // type == -1 : query the min of all the types of diagrams.
683 // type = 0 : min, 1 : sad, 2 : max
684
685 double min_persistence = std::numeric_limits<double>::max();
686 if(do_min_ && (type == -1 || type == 0)) {
687 for(size_t i = 0; i < bidder_diagrams_min_.size(); ++i) {
688 for(size_t j = 0; j < bidder_diagrams_min_[i].size(); ++j) {
689 Bidder b = bidder_diagrams_min_[i].at(j);
690 double persistence = b.getPersistence();
691 if(persistence < min_persistence) {
692 min_persistence = persistence;
693 }
694 }
695 }
696 }
697
698 if(do_sad_ && (type == -1 || type == 1)) {
699 for(size_t i = 0; i < bidder_diagrams_saddle_.size(); ++i) {
700 for(size_t j = 0; j < bidder_diagrams_saddle_[i].size(); ++j) {
701 Bidder b = bidder_diagrams_saddle_[i].at(j);
702 double persistence = b.getPersistence();
703 if(persistence < min_persistence) {
704 min_persistence = persistence;
705 }
706 }
707 }
708 }
709
710 if(do_max_ && (type == -1 || type == 2)) {
711 for(size_t i = 0; i < bidder_diagrams_max_.size(); ++i) {
712 for(size_t j = 0; j < bidder_diagrams_max_[i].size(); ++j) {
713 Bidder b = bidder_diagrams_max_[i].at(j);
714 double persistence = b.getPersistence();
715 if(persistence < min_persistence) {
716 min_persistence = persistence;
717 }
718 }
719 }
720 }
721 return min_persistence;
722}
723
724std::vector<std::vector<double>> ttk::PDClustering::getMinPrices() {
725 std::vector<std::vector<double>> min_prices(3);
726
727 if(original_dos[0]) {
728 for(int i = 0; i < numberOfInputs_; ++i) {
729 min_prices[0].emplace_back(std::numeric_limits<double>::max());
730 for(size_t j = 0; j < centroids_with_price_min_[i].size(); ++j) {
731 Good g = centroids_with_price_min_[i].at(j);
732 double price = g.getPrice();
733 if(price < min_prices[0][i]) {
734 min_prices[0][i] = price;
735 }
736 }
737 }
738 }
739
740 if(original_dos[1]) {
741 for(int i = 0; i < numberOfInputs_; ++i) {
742 min_prices[1].emplace_back(std::numeric_limits<double>::max());
743 for(size_t j = 0; j < centroids_with_price_saddle_[i].size(); ++j) {
744 Good g = centroids_with_price_saddle_[i].at(j);
745 double price = g.getPrice();
746 if(price < min_prices[1][i]) {
747 min_prices[1][i] = price;
748 }
749 }
750 }
751 }
752
753 if(original_dos[2]) {
754 for(int i = 0; i < numberOfInputs_; ++i) {
755 min_prices[2].emplace_back(std::numeric_limits<double>::max());
756 for(size_t j = 0; j < centroids_with_price_max_[i].size(); ++j) {
757 Good g = centroids_with_price_max_[i].at(j);
758 double price = g.getPrice();
759 if(price < min_prices[2][i]) {
760 min_prices[2][i] = price;
761 }
762 }
763 }
764 }
765
766 return min_prices;
767}
768
769std::vector<std::vector<double>> ttk::PDClustering::getMinDiagonalPrices() {
770 std::vector<std::vector<double>> min_prices(3);
771 if(original_dos[0]) {
772 for(int i = 0; i < numberOfInputs_; ++i) {
773 min_prices[0].emplace_back(std::numeric_limits<double>::max());
774 for(size_t j = 0; j < current_bidder_diagrams_min_[i].size(); ++j) {
775 Bidder b = current_bidder_diagrams_min_[i].at(j);
776 double price = b.diagonal_price_;
777 if(price < min_prices[0][i]) {
778 min_prices[0][i] = price;
779 }
780 }
781 if(min_prices[0][i] >= std::numeric_limits<double>::max() / 2.) {
782 min_prices[0][i] = 0;
783 }
784 }
785 }
786
787 if(original_dos[1]) {
788 for(int i = 0; i < numberOfInputs_; ++i) {
789 min_prices[1].emplace_back(std::numeric_limits<double>::max());
790 for(size_t j = 0; j < current_bidder_diagrams_saddle_[i].size(); ++j) {
791 Bidder b = current_bidder_diagrams_saddle_[i].at(j);
792 double price = b.diagonal_price_;
793 if(price < min_prices[1][i]) {
794 min_prices[1][i] = price;
795 }
796 }
797 if(min_prices[1][i] >= std::numeric_limits<double>::max() / 2.) {
798 min_prices[1][i] = 0;
799 }
800 }
801 }
802
803 if(original_dos[2]) {
804 for(int i = 0; i < numberOfInputs_; ++i) {
805 min_prices[2].emplace_back(std::numeric_limits<double>::max());
806 for(size_t j = 0; j < current_bidder_diagrams_max_[i].size(); ++j) {
807 Bidder b = current_bidder_diagrams_max_[i].at(j);
808 double price = b.diagonal_price_;
809 if(price < min_prices[2][i]) {
810 min_prices[2][i] = price;
811 }
812 }
813 if(min_prices[2][i] >= std::numeric_limits<double>::max() / 2.) {
814 min_prices[2][i] = 0;
815 }
816 }
817 }
818 return min_prices;
819}
820
822 const BidderDiagram &D2,
823 const double delta_lim) {
824 GoodDiagram D2_bis = diagramToCentroid(D2);
825 return computeDistance(D1, D2_bis, delta_lim);
826}
827
829 const GoodDiagram &D2,
830 const double delta_lim) {
831 std::vector<MatchingType> matchings;
832 const auto D2_bis = centroidWithZeroPrices(D2);
834 wasserstein_, geometrical_factor_, lambda_, delta_lim, use_kdtree_);
835 auction.BuildAuctionDiagrams(D1, D2_bis);
836 double cost = auction.run(matchings);
837 return cost;
838}
839
841 const GoodDiagram *const D2,
842 const double delta_lim) {
843 std::vector<MatchingType> matchings;
845 wasserstein_, geometrical_factor_, lambda_, delta_lim, use_kdtree_);
846 int size1 = D1->size();
847 auction.BuildAuctionDiagrams(*D1, *D2);
848 double cost = auction.run(matchings);
849 // Diagonal Points were added in the original diagram. The following line
850 // removes them.
851 D1->resize(size1);
852 return cost;
853}
854
856 const GoodDiagram &D2,
857 const double delta_lim) {
858 BidderDiagram D1_bis = centroidToDiagram(D1);
859 return computeDistance(D1_bis, D2, delta_lim);
860}
861
865 for(size_t i = 0; i < centroid.size(); i++) {
866 Good g = centroid.at(i);
867 g.setPrice(0);
868 GD.emplace_back(g);
869 }
870 return GD;
871}
872
876 for(size_t i = 0; i < diagram.size(); i++) {
877 Bidder b = diagram.at(i);
878 b.setDiagonalPrice(0);
879 BD.emplace_back(b);
880 }
881 return BD;
882}
883
887 for(size_t i = 0; i < centroid.size(); i++) {
888 Good g = centroid.at(i);
889
890 Bidder b = Bidder(g.x_, g.y_, g.isDiagonal(), BD.size());
891 b.SetCriticalCoordinates(g.coords_[0], g.coords_[1], g.coords_[2]);
892 b.setPositionInAuction(BD.size());
893 BD.emplace_back(b);
894 }
895 return BD;
896}
897
901 for(size_t i = 0; i < diagram.size(); i++) {
902 Bidder b = diagram.at(i);
903
904 Good g = Good(b.x_, b.y_, b.isDiagonal(), GD.size());
905 g.SetCriticalCoordinates(b.coords_[0], b.coords_[1], b.coords_[2]);
906 GD.emplace_back(g);
907 }
908 return GD;
909}
910
912 clustering_ = std::vector<std::vector<int>>(k_);
913}
914
916 std::vector<int> idx(numberOfInputs_);
917 // To perform a random draw with replacement, the vector {1, 2, ...,
918 // numberOfInputs_} is shuffled, and we consider its k_ first elements to be
919 // the initial centroids.
920 for(int i = 0; i < numberOfInputs_; i++) {
921 idx[i] = i;
922 }
923 if(!deterministic_)
924 std::shuffle(idx.begin(), idx.end(), std::random_device());
925
926 for(int c = 0; c < k_; c++) {
927 if(do_min_) {
928 GoodDiagram centroid_min
929 = diagramToCentroid(current_bidder_diagrams_min_[idx[c]]);
930 centroids_min_.emplace_back(centroid_min);
931 }
932 if(do_sad_) {
933 GoodDiagram centroid_sad
934 = diagramToCentroid(current_bidder_diagrams_saddle_[idx[c]]);
935 centroids_saddle_.emplace_back(centroid_sad);
936 }
937 if(do_max_) {
938 GoodDiagram centroid_max
939 = diagramToCentroid(current_bidder_diagrams_max_[idx[c]]);
940 centroids_max_.emplace_back(centroid_max);
941 }
942 }
943}
944
946 std::vector<int> indexes_clusters;
947 int random_idx = deterministic_ ? 0 : rand() % numberOfInputs_;
948 indexes_clusters.emplace_back(random_idx);
949
950 if(do_min_) {
951 GoodDiagram centroid_min
952 = diagramToCentroid(current_bidder_diagrams_min_[random_idx]);
953 centroids_min_.emplace_back(centroid_min);
954 }
955 if(do_sad_) {
956 GoodDiagram centroid_sad
957 = diagramToCentroid(current_bidder_diagrams_saddle_[random_idx]);
958 centroids_saddle_.emplace_back(centroid_sad);
959 }
960 if(do_max_) {
961 GoodDiagram centroid_max
962 = diagramToCentroid(current_bidder_diagrams_max_[random_idx]);
963 centroids_max_.emplace_back(centroid_max);
964 }
965
966 const auto square = [](const double a) { return a * a; };
967
968 while((int)indexes_clusters.size() < k_) {
969 std::vector<double> min_distance_to_centroid(numberOfInputs_);
970 std::vector<double> probabilities(numberOfInputs_);
971
972 // Uncomment for a deterministic algorithm
973 double maximal_distance = 0;
974 int candidate_centroid = 0;
975
976 for(int i = 0; i < numberOfInputs_; i++) {
977
978 min_distance_to_centroid[i] = std::numeric_limits<double>::max();
979 if(std::find(indexes_clusters.begin(), indexes_clusters.end(), i)
980 != indexes_clusters.end()) {
981
982 min_distance_to_centroid[i] = 0;
983 } else {
984
985 for(size_t j = 0; j < indexes_clusters.size(); ++j) {
986
987 double distance = 0;
988 if(do_min_) {
989 GoodDiagram centroid_min
990 = centroidWithZeroPrices(centroids_min_[j]);
991 distance += computeDistance(
992 current_bidder_diagrams_min_[i], centroid_min, 0.01);
993 }
994 if(do_sad_) {
995 GoodDiagram centroid_saddle
996 = centroidWithZeroPrices(centroids_saddle_[j]);
997 distance += computeDistance(
998 current_bidder_diagrams_saddle_[i], centroid_saddle, 0.01);
999 }
1000 if(do_max_) {
1001 GoodDiagram centroid_max
1002 = centroidWithZeroPrices(centroids_max_[j]);
1003 distance += computeDistance(
1004 current_bidder_diagrams_max_[i], centroid_max, 0.01);
1005 }
1006 if(distance < min_distance_to_centroid[i]) {
1007 min_distance_to_centroid[i] = distance;
1008 }
1009 }
1010 }
1011 probabilities[i] = square(min_distance_to_centroid[i]);
1012
1013 // The following block is useful in case of need for a deterministic
1014 // algorithm
1015 if(deterministic_ && min_distance_to_centroid[i] > maximal_distance) {
1016 maximal_distance = min_distance_to_centroid[i];
1017 candidate_centroid = i;
1018 }
1019 }
1020 // Comment the following four lines to make it deterministic
1021 std::random_device rd;
1022 std::mt19937 gen(rd());
1023 std::discrete_distribution<int> distribution(
1024 probabilities.begin(), probabilities.end());
1025
1026 if(!deterministic_) {
1027 candidate_centroid = distribution(gen);
1028 }
1029
1030 indexes_clusters.emplace_back(candidate_centroid);
1031 if(do_min_) {
1032 GoodDiagram centroid_min
1033 = diagramToCentroid(current_bidder_diagrams_min_[candidate_centroid]);
1034 centroids_min_.emplace_back(centroid_min);
1035 }
1036 if(do_sad_) {
1037 GoodDiagram centroid_sad = diagramToCentroid(
1038 current_bidder_diagrams_saddle_[candidate_centroid]);
1039 centroids_saddle_.emplace_back(centroid_sad);
1040 }
1041 if(do_max_) {
1042 GoodDiagram centroid_max
1043 = diagramToCentroid(current_bidder_diagrams_max_[candidate_centroid]);
1044 centroids_max_.emplace_back(centroid_max);
1045 }
1046 }
1047}
1048
1050 // r_ is a vector stating for each diagram if its distance to its centroid is
1051 // up to date (false) or needs to be recomputed (true)
1052 r_ = std::vector<bool>(numberOfInputs_);
1053 // u_ is a vector of upper bounds of the distance of each diagram to its
1054 // closest centroid
1055 u_ = std::vector<double>(numberOfInputs_);
1056 inv_clustering_ = std::vector<int>(numberOfInputs_);
1057 for(int i = 0; i < numberOfInputs_; i++) {
1058 r_[i] = true;
1059 u_[i] = std::numeric_limits<double>::max();
1060 inv_clustering_[i] = -1;
1061 }
1062 // l_ is the matrix of lower bounds for the distance from each diagram
1063 // to each centroid
1064 l_ = std::vector<std::vector<double>>(numberOfInputs_);
1065 for(int i = 0; i < numberOfInputs_; ++i) {
1066 l_[i] = std::vector<double>(k_);
1067 for(int c = 0; c < k_; ++c) {
1068 l_[i][c] = 0;
1069 }
1070 }
1071
1072 // And d_ is a (K x K) matrix storing the distances between each pair of
1073 // centroids
1074 centroidsDistanceMatrix_.resize(k_);
1075 for(int i = 0; i < k_; ++i) {
1076 centroidsDistanceMatrix_[i].resize(k_, 0.0);
1077 }
1078}
1079
1080std::vector<std::vector<double>> ttk::PDClustering::getDistanceMatrix() {
1081 std::vector<std::vector<double>> D(numberOfInputs_);
1082
1083 for(int i = 0; i < numberOfInputs_; ++i) {
1084 BidderDiagram D1_min, D1_sad, D1_max;
1085 if(do_min_) {
1086 D1_min = diagramWithZeroPrices(current_bidder_diagrams_min_[i]);
1087 }
1088 if(do_sad_) {
1089 D1_sad = diagramWithZeroPrices(current_bidder_diagrams_saddle_[i]);
1090 }
1091 if(do_max_) {
1092 D1_max = diagramWithZeroPrices(current_bidder_diagrams_max_[i]);
1093 }
1094 for(int c = 0; c < k_; ++c) {
1095 GoodDiagram D2_min, D2_sad, D2_max;
1096 double distance = 0;
1097 if(do_min_) {
1098 D2_min = centroids_min_[c];
1099 distance += computeDistance(D1_min, D2_min, 0.01);
1100 }
1101 if(do_sad_) {
1102 D2_sad = centroids_saddle_[c];
1103 distance += computeDistance(D1_sad, D2_sad, 0.01);
1104 }
1105 if(do_max_) {
1106 D2_max = centroids_max_[c];
1107 distance += computeDistance(D1_max, D2_max, 0.01);
1108 }
1109 D[i].emplace_back(distance);
1110 }
1111 }
1112 return D;
1113}
1114
1116 for(int i = 0; i < k_; ++i) {
1117 GoodDiagram D1_min, D1_sad, D1_max;
1118 if(do_min_) {
1119 D1_min = centroidWithZeroPrices(centroids_min_[i]);
1120 }
1121 if(do_sad_) {
1122 D1_sad = centroidWithZeroPrices(centroids_saddle_[i]);
1123 }
1124 if(do_max_) {
1125 D1_max = centroidWithZeroPrices(centroids_max_[i]);
1126 }
1127 for(int j = i + 1; j < k_; ++j) {
1128 double distance{};
1129 GoodDiagram D2_min, D2_sad, D2_max;
1130 if(do_min_) {
1131 D2_min = centroidWithZeroPrices(centroids_min_[j]);
1132 distance += computeDistance(D1_min, D2_min, 0.01);
1133 }
1134 if(do_sad_) {
1135 D2_sad = centroidWithZeroPrices(centroids_saddle_[j]);
1136 distance += computeDistance(D1_sad, D2_sad, 0.01);
1137 }
1138 if(do_max_) {
1139 D2_max = centroidWithZeroPrices(centroids_max_[j]);
1140 distance += computeDistance(D1_max, D2_max, 0.01);
1141 }
1142
1143 centroidsDistanceMatrix_[i][j] = distance;
1144 centroidsDistanceMatrix_[j][i] = distance;
1145 }
1146 }
1147}
1148
1150 distanceToCentroid_.resize(numberOfInputs_);
1151
1152 for(int i = 0; i < numberOfInputs_; ++i) {
1153 double delta_lim{0.01};
1154 double distance{};
1155 auto c = inv_clustering_[i];
1156 if(original_dos[0]) {
1157 GoodDiagram centroid_min = centroidWithZeroPrices(centroids_min_[c]);
1158 BidderDiagram bidder_diag
1159 = diagramWithZeroPrices(current_bidder_diagrams_min_[i]);
1160 distance += computeDistance(bidder_diag, centroid_min, delta_lim);
1161 }
1162 if(original_dos[1]) {
1163 GoodDiagram centroid_saddle
1164 = centroidWithZeroPrices(centroids_saddle_[c]);
1165 BidderDiagram bidder_diag
1166 = diagramWithZeroPrices(current_bidder_diagrams_saddle_[i]);
1167 distance += computeDistance(bidder_diag, centroid_saddle, delta_lim);
1168 }
1169 if(original_dos[2]) {
1170 GoodDiagram centroid_max = centroidWithZeroPrices(centroids_max_[c]);
1171 BidderDiagram bidder_diag
1172 = diagramWithZeroPrices(current_bidder_diagrams_max_[i]);
1173 distance += computeDistance(bidder_diag, centroid_max, delta_lim);
1174 }
1175 distanceToCentroid_[i] = distance;
1176 }
1177}
1178
1180 if(k_ > 1) {
1181 std::vector<std::vector<double>> distance_matrix = getDistanceMatrix();
1182 old_clustering_ = clustering_;
1183 invertClusters();
1184 initializeEmptyClusters();
1185
1186 for(int i = 0; i < numberOfInputs_; ++i) {
1187 double min_distance_to_centroid = std::numeric_limits<double>::max();
1188 int cluster = -1;
1189 for(int c = 0; c < k_; ++c) {
1190 if(distance_matrix[i][c] < min_distance_to_centroid) {
1191 min_distance_to_centroid = distance_matrix[i][c];
1192 cluster = c;
1193 }
1194 }
1195
1196 clustering_[cluster].emplace_back(i);
1197 if(cluster != inv_clustering_[i]) {
1198 // New centroid attributed to this diagram
1199 resetDosToOriginalValues();
1200 barycenter_inputs_reset_flag = true;
1201 if(do_min_) {
1202 centroids_with_price_min_[i]
1203 = centroidWithZeroPrices(centroids_min_[cluster]);
1204 }
1205 if(do_sad_) {
1206 centroids_with_price_saddle_[i]
1207 = centroidWithZeroPrices(centroids_saddle_[cluster]);
1208 }
1209 if(do_max_) {
1210 centroids_with_price_max_[i]
1211 = centroidWithZeroPrices(centroids_max_[cluster]);
1212 }
1213 inv_clustering_[i] = cluster;
1214 }
1215 }
1216
1217 } else {
1218 old_clustering_ = clustering_;
1219 invertClusters();
1220 initializeEmptyClusters();
1221
1222 for(int i = 0; i < numberOfInputs_; i++) {
1223 clustering_[0].emplace_back(i);
1224 if(n_iterations_ < 1) {
1225 if(do_min_) {
1226 centroids_with_price_min_[i]
1227 = centroidWithZeroPrices(centroids_min_[0]);
1228 }
1229 if(do_sad_) {
1230 centroids_with_price_saddle_[i]
1231 = centroidWithZeroPrices(centroids_saddle_[0]);
1232 }
1233 if(do_max_) {
1234 centroids_with_price_max_[i]
1235 = centroidWithZeroPrices(centroids_max_[0]);
1236 }
1237 }
1238 inv_clustering_[i] = 0;
1239 }
1240 }
1241}
1242
1247
1248 // Initializes clusters with -1
1249 inv_clustering_ = std::vector<int>(numberOfInputs_);
1250 for(int i = 0; i < numberOfInputs_; ++i) {
1251 inv_clustering_[i] = -1;
1252 }
1253
1254 // Fill in the clusters
1255 for(int c = 0; c < k_; ++c) {
1256 for(size_t j = 0; j < clustering_[c].size(); ++j) {
1257 int idx = clustering_[c][j];
1258 inv_clustering_[idx] = c;
1259 }
1260 }
1261}
1262
1264 clustering_ = std::vector<std::vector<int>>(k_);
1265 for(int i = 0; i < numberOfInputs_; ++i) {
1266 clustering_[inv_clustering_[i]].emplace_back(i);
1267 }
1268
1269 // Check if a cluster was left without diagram
1270 for(int c = 0; c < k_; ++c) {
1271 if(clustering_[c].empty()) {
1272 std::cout << "Problem in invertInverseClusters()... \nCluster " << c
1273 << " was left with no diagram attached to it... " << std::endl;
1274 }
1275 }
1276}
1277
1279 // Step 1
1280 getCentroidDistanceMatrix();
1281 old_clustering_ = clustering_;
1282 // self.old_clusters = copy.copy(self.clusters)
1283 invertClusters();
1284 initializeEmptyClusters();
1285 bool do_min = original_dos[0];
1286 bool do_sad = original_dos[1];
1287 bool do_max = original_dos[2];
1288
1289 for(int i = 0; i < numberOfInputs_; ++i) {
1290 // Step 3 find potential changes of clusters
1291 BidderDiagram D1_min, D1_sad, D1_max;
1292 if(do_min) {
1293 D1_min = diagramWithZeroPrices(current_bidder_diagrams_min_[i]);
1294 }
1295 if(do_sad) {
1296 D1_sad = diagramWithZeroPrices(current_bidder_diagrams_saddle_[i]);
1297 }
1298 if(do_max) {
1299 D1_max = diagramWithZeroPrices(current_bidder_diagrams_max_[i]);
1300 }
1301
1302 for(int c = 0; c < k_; ++c) {
1303 if(inv_clustering_[i] == -1) {
1304 // If not yet assigned, assign it first to a random cluster
1305
1306 if(deterministic_) {
1307 inv_clustering_[i] = i % k_;
1308 } else {
1309 std::cout << " - ASSIGNED TO A RANDOM CLUSTER " << '\n';
1310 inv_clustering_[i] = rand() % (k_);
1311 }
1312
1313 r_[i] = true;
1314 if(do_min) {
1315 centroids_with_price_min_[i]
1316 = centroidWithZeroPrices(centroids_min_[inv_clustering_[i]]);
1317 }
1318 if(do_sad) {
1319 centroids_with_price_saddle_[i]
1320 = centroidWithZeroPrices(centroids_saddle_[inv_clustering_[i]]);
1321 }
1322 if(do_max) {
1323 centroids_with_price_max_[i]
1324 = centroidWithZeroPrices(centroids_max_[inv_clustering_[i]]);
1325 }
1326 }
1327
1328 if(c != inv_clustering_[i] && u_[i] > l_[i][c]
1329 && u_[i] > 0.5 * centroidsDistanceMatrix_[inv_clustering_[i]][c]) {
1330 // Step 3a, If necessary, recompute the distance to centroid
1331 if(r_[i]) {
1332 double distance = 0;
1333 GoodDiagram centroid_min, centroid_sad, centroid_max;
1334 if(do_min) {
1335 centroid_min
1336 = centroidWithZeroPrices(centroids_min_[inv_clustering_[i]]);
1337 distance += computeDistance(D1_min, centroid_min, 0.01);
1338 }
1339 if(do_sad) {
1340 centroid_sad
1341 = centroidWithZeroPrices(centroids_saddle_[inv_clustering_[i]]);
1342 distance += computeDistance(D1_sad, centroid_sad, 0.01);
1343 }
1344 if(do_max) {
1345 centroid_max
1346 = centroidWithZeroPrices(centroids_max_[inv_clustering_[i]]);
1347 distance += computeDistance(D1_max, centroid_max, 0.01);
1348 }
1349 r_[i] = false;
1350 u_[i] = distance;
1351 l_[i][inv_clustering_[i]] = distance;
1352 }
1353 // Step 3b, check if still potential change of clusters
1354 if((n_iterations_ > 2 || n_iterations_ < 1)
1355 && (u_[i] > l_[i][c]
1356 || u_[i]
1357 > 0.5 * centroidsDistanceMatrix_[inv_clustering_[i]][c])) {
1358 BidderDiagram diagram_min, diagram_sad, diagram_max;
1359 GoodDiagram centroid_min, centroid_sad, centroid_max;
1360 double distance = 0;
1361
1362 if(do_min) {
1363 centroid_min = centroidWithZeroPrices(centroids_min_[c]);
1364 diagram_min
1365 = diagramWithZeroPrices(current_bidder_diagrams_min_[i]);
1366 distance += computeDistance(diagram_min, centroid_min, 0.01);
1367 }
1368 if(do_sad) {
1369 centroid_sad = centroidWithZeroPrices(centroids_saddle_[c]);
1370 diagram_sad
1371 = diagramWithZeroPrices(current_bidder_diagrams_saddle_[i]);
1372 distance += computeDistance(diagram_sad, centroid_sad, 0.01);
1373 }
1374 if(do_max) {
1375 centroid_max = centroidWithZeroPrices(centroids_max_[c]);
1376 diagram_max
1377 = diagramWithZeroPrices(current_bidder_diagrams_max_[i]);
1378 distance += computeDistance(diagram_max, centroid_max, 0.01);
1379 }
1380 l_[i][c] = distance;
1381 // TODO Prices are lost here... If distance<self.u[i], we should keep
1382 // the prices
1383 if(distance < u_[i]) {
1384 // Changing cluster
1385 resetDosToOriginalValues();
1386 barycenter_inputs_reset_flag = true;
1387 u_[i] = distance;
1388 inv_clustering_[i] = c;
1389
1390 if(do_min) {
1391 centroids_with_price_min_[i]
1392 = centroidWithZeroPrices(centroids_min_[c]);
1393 }
1394 if(do_sad) {
1395 centroids_with_price_saddle_[i]
1396 = centroidWithZeroPrices(centroids_saddle_[c]);
1397 }
1398 if(do_max) {
1399 centroids_with_price_max_[i]
1400 = centroidWithZeroPrices(centroids_max_[c]);
1401 }
1402 }
1403 }
1404 }
1405 }
1406 }
1407 invertInverseClusters();
1408 for(int c = 0; c < k_; ++c) {
1409 if(clustering_[c].empty()) {
1410 std::cout << "Adding artificial centroid because a cluster was empty"
1411 << std::endl;
1412 bool idx_acceptable = false;
1413 int idx = 0;
1414
1415 std::vector<double> copy_of_u(u_.size());
1416 copy_of_u = u_;
1417 while(!idx_acceptable) {
1418 auto argMax = std::max_element(copy_of_u.begin(), copy_of_u.end());
1419 idx = std::distance(copy_of_u.begin(), argMax);
1420 if(inv_clustering_[idx] < k_ && inv_clustering_[idx] >= 0
1421 && clustering_[inv_clustering_[idx]].size() > 1) {
1422 idx_acceptable = true;
1423 int cluster_removal = inv_clustering_[idx];
1424 // Removing the index to remove
1425 clustering_[cluster_removal].erase(
1426 std::remove(clustering_[cluster_removal].begin(),
1427 clustering_[cluster_removal].end(), idx),
1428 clustering_[cluster_removal].end());
1429 } else {
1430 if(copy_of_u.size() > (size_t)idx) {
1431 copy_of_u.erase(argMax);
1432 } else {
1433 idx_acceptable = true;
1434 int cluster_max = 0;
1435 if(!clustering_[cluster_max].empty()) {
1436 idx = clustering_[cluster_max][0];
1437 }
1438 for(int i_test = 1; i_test < k_; i_test++) {
1439 if(clustering_[i_test].size() > clustering_[cluster_max].size()) {
1440 cluster_max = i_test;
1441 idx = clustering_[cluster_max][0];
1442 }
1443 }
1444 int cluster_removal = inv_clustering_[idx];
1445 clustering_[cluster_removal].erase(
1446 std::remove(clustering_[cluster_removal].begin(),
1447 clustering_[cluster_removal].end(), idx),
1448 clustering_[cluster_removal].end());
1449 }
1450 }
1451 }
1452
1453 clustering_[c].emplace_back(idx);
1454 inv_clustering_[idx] = c;
1455
1456 if(do_min) {
1457 centroids_min_[c]
1458 = diagramToCentroid(current_bidder_diagrams_min_[idx]);
1459 centroids_with_price_min_[idx]
1460 = centroidWithZeroPrices(centroids_min_[c]);
1461 }
1462 if(do_sad) {
1463 centroids_saddle_[c]
1464 = diagramToCentroid(current_bidder_diagrams_saddle_[idx]);
1465 centroids_with_price_saddle_[idx]
1466 = centroidWithZeroPrices(centroids_saddle_[c]);
1467 }
1468 if(do_max) {
1469 centroids_max_[c]
1470 = diagramToCentroid(current_bidder_diagrams_max_[idx]);
1471 centroids_with_price_max_[idx]
1472 = centroidWithZeroPrices(centroids_max_[c]);
1473 }
1474 resetDosToOriginalValues();
1475 barycenter_inputs_reset_flag = true;
1476 }
1477 }
1478}
1479
1481 std::vector<std::vector<double>> *min_price,
1482 std::vector<std::vector<double>> *min_diag_price,
1483 std::vector<std::vector<std::vector<std::vector<MatchingType>>>>
1484 &all_matchings_per_type_and_cluster,
1485 int only_matchings) {
1486 barycenter_inputs_reset_flag = true;
1487 std::vector<double> max_shift_vector(3);
1488 max_shift_vector[0] = 0;
1489 max_shift_vector[1] = 0;
1490 max_shift_vector[2] = 0;
1491 double max_shift_c_min = 0;
1492 double max_shift_c_sad = 0;
1493 double max_shift_c_max = 0;
1494 double max_wasserstein_shift = 0;
1495 bool precision_min = true;
1496 bool precision_sad = true;
1497 bool precision_max = true;
1498 cost_ = 0;
1499 double sq_dist_min = cost_min_;
1500 double sq_dist_sad = cost_sad_;
1501 double sq_dist_max = cost_max_;
1502 if(do_min_) {
1503 cost_min_ = 0;
1504 }
1505 if(do_sad_) {
1506 cost_sad_ = 0;
1507 }
1508 if(do_max_) {
1509 cost_max_ = 0;
1510 }
1511
1512 for(int c = 0; c < k_; ++c) {
1513 if(!clustering_[c].empty()) {
1514 std::vector<GoodDiagram> centroids_with_price_min,
1515 centroids_with_price_sad, centroids_with_price_max;
1516 int count = 0;
1517 for(int idx : clustering_[c]) {
1518 // Timer time_first_thing;
1519 // Find the position of diagrams[idx] in old cluster c
1520 std::vector<int>::iterator i = std::find(
1521 old_clustering_[c].begin(), old_clustering_[c].end(), idx);
1522 int pos = (i == old_clustering_[c].end())
1523 ? -1
1524 : std::distance(old_clustering_[c].begin(), i);
1525 if(pos >= 0) {
1526 // Diagram was already linked to this centroid before
1527 if(do_min_) {
1528 centroids_with_price_min.emplace_back(
1529 centroids_with_price_min_[idx]);
1530 }
1531 if(do_sad_) {
1532 centroids_with_price_sad.emplace_back(
1533 centroids_with_price_saddle_[idx]);
1534 }
1535 if(do_max_) {
1536 centroids_with_price_max.emplace_back(
1537 centroids_with_price_max_[idx]);
1538 }
1539 } else {
1540 int number_of_points_min = 0;
1541 int number_of_points_max = 0;
1542 int number_of_points_sad = 0;
1543
1544 // Otherwise, centroid is given 0 prices and the diagram is given 0
1545 // diagonal-prices
1546 if(do_min_) {
1547 centroids_with_price_min.emplace_back(
1548 centroidWithZeroPrices(centroids_min_[c]));
1549 current_bidder_diagrams_min_[idx]
1550 = diagramWithZeroPrices(current_bidder_diagrams_min_[idx]);
1551 number_of_points_min += centroids_with_price_min_[idx].size()
1552 + current_bidder_diagrams_min_[idx].size();
1553 }
1554 if(do_sad_) {
1555 centroids_with_price_sad.emplace_back(
1556 centroidWithZeroPrices(centroids_saddle_[c]));
1557 current_bidder_diagrams_saddle_[idx]
1558 = diagramWithZeroPrices(current_bidder_diagrams_saddle_[idx]);
1559 number_of_points_sad
1560 += centroids_with_price_saddle_[idx].size()
1561 + current_bidder_diagrams_saddle_[idx].size();
1562 }
1563 if(do_max_) {
1564 centroids_with_price_max.emplace_back(
1565 centroidWithZeroPrices(centroids_max_[c]));
1566 current_bidder_diagrams_max_[idx]
1567 = diagramWithZeroPrices(current_bidder_diagrams_max_[idx]);
1568 number_of_points_max += centroids_with_price_max_[idx].size()
1569 + current_bidder_diagrams_max_[idx].size();
1570 }
1571
1572 if(n_iterations_ > 1) {
1573 // If diagram new to cluster and we're not at first iteration,
1574 // precompute prices for the objects via compute_distance()
1575 // number_of_points /= (int)do_min_ + (int)do_sad_ + (int)do_max_;
1576 // double d_estimated = pow(cost_ / numberOfInputs_, 1. /
1577 // wasserstein_) + 1e-7; We use pointer in the auction in order to
1578 // keep the prices at the end
1579
1580 if(do_min_) {
1581 // double estimated_delta_lim = number_of_points_min *
1582 // epsilon_[0] / (2*sq_dist_min) ;
1583 double estimated_delta_lim
1584 = 1.
1585 / sqrt(1 - number_of_points_min * epsilon_[0] / sq_dist_min)
1586 - 1;
1587 if(estimated_delta_lim > 1) {
1588 estimated_delta_lim = 1;
1589 }
1590 computeDistance(&(current_bidder_diagrams_min_[idx]),
1591 &(centroids_with_price_min[count]),
1592 estimated_delta_lim);
1593 }
1594 if(do_sad_) {
1595 // double estimated_delta_lim = number_of_points_sad *
1596 // epsilon_[1] / (2*sq_dist_sad);
1597 double estimated_delta_lim
1598 = 1.
1599 / sqrt(1 - number_of_points_sad * epsilon_[1] / sq_dist_sad)
1600 - 1;
1601 if(estimated_delta_lim > 1) {
1602 estimated_delta_lim = 1;
1603 }
1604 computeDistance(&(current_bidder_diagrams_saddle_[idx]),
1605 &(centroids_with_price_sad[count]),
1606 estimated_delta_lim);
1607 }
1608 if(do_max_) {
1609 // double estimated_delta_lim = number_of_points_max *
1610 // epsilon_[2] / (2*sq_dist_max);
1611 double estimated_delta_lim
1612 = 1.
1613 / sqrt(1 - number_of_points_max * epsilon_[2] / sq_dist_max)
1614 - 1;
1615 if(estimated_delta_lim > 1) {
1616 estimated_delta_lim = 1;
1617 }
1618 computeDistance(&(current_bidder_diagrams_max_[idx]),
1619 &(centroids_with_price_max[count]),
1620 estimated_delta_lim);
1621 }
1622 }
1623 }
1624 count++;
1625 }
1626 double total_cost = 0;
1627 double wasserstein_shift = 0;
1628
1629 using KDTreePair = PDBarycenter::KDTreePair;
1630
1631 if(do_min_) {
1632 std::vector<std::vector<MatchingType>> all_matchings;
1633 std::vector<int> sizes;
1634 Timer time_preprocess_bary;
1635 std::vector<BidderDiagram> diagrams_c_min;
1636 if(barycenter_inputs_reset_flag) {
1637 for(int idx : clustering_[c]) {
1638 diagrams_c_min.emplace_back(current_bidder_diagrams_min_[idx]);
1639 }
1640 sizes.resize(diagrams_c_min.size());
1641 for(size_t i = 0; i < diagrams_c_min.size(); i++) {
1642 sizes[i] = diagrams_c_min[i].size();
1643 }
1644 barycenter_computer_min_[c].setNumberOfInputs(diagrams_c_min.size());
1645 barycenter_computer_min_[c].setCurrentBidders(diagrams_c_min);
1646
1647 std::vector<GoodDiagram> barycenter_goods(clustering_[c].size());
1648 for(size_t i_diagram = 0; i_diagram < clustering_[c].size();
1649 i_diagram++) {
1650 barycenter_goods[i_diagram]
1651 = centroids_with_price_min_[clustering_[c][i_diagram]];
1652 }
1653 barycenter_computer_min_[c].setCurrentBarycenter(barycenter_goods);
1654 all_matchings.resize(diagrams_c_min.size());
1655 } else {
1656 sizes.resize(barycenter_computer_min_[c].getCurrentBidders().size());
1657 for(size_t i = 0;
1658 i < barycenter_computer_min_[c].getCurrentBidders().size(); i++) {
1659 sizes[i]
1660 = barycenter_computer_min_[c].getCurrentBidders().at(i).size();
1661 }
1662 diagrams_c_min.resize(
1663 barycenter_computer_min_[c].getCurrentBidders().size());
1664 all_matchings.resize(
1665 barycenter_computer_min_[c].getCurrentBidders().size());
1666 }
1667
1668 KDTreePair pair;
1669 bool use_kdt = false;
1670 if(!barycenter_computer_min_[c].getCurrentBarycenter()[0].empty()) {
1671 pair = barycenter_computer_min_[c].getKDTree();
1672 use_kdt = true;
1673 }
1674
1675 barycenter_computer_min_[c].runMatching(
1676 &total_cost, epsilon_[0], sizes, *pair.first, pair.second,
1677 &(min_diag_price->at(0)), &(min_price->at(0)), &(all_matchings),
1678 use_kdt, only_matchings);
1679 for(size_t ii = 0; ii < all_matchings.size(); ii++) {
1680 all_matchings_per_type_and_cluster[c][0][ii].resize(
1681 all_matchings[ii].size());
1682 all_matchings_per_type_and_cluster[c][0][ii] = all_matchings[ii];
1683 }
1684 for(int ii = all_matchings.size(); ii < numberOfInputs_; ii++) {
1685 all_matchings_per_type_and_cluster[c][0][ii].resize(0);
1686 }
1687
1688 precision_min
1689 = barycenter_computer_min_[c].isPrecisionObjectiveMet(deltaLim_, 0);
1690 cost_min_ += total_cost;
1691 Timer time_update;
1692 if(!only_matchings) {
1693 max_shift_c_min
1694 = barycenter_computer_min_[c].updateBarycenter(all_matchings);
1695 }
1696
1697 if(max_shift_c_min > max_shift_vector[0]) {
1698 max_shift_vector[0] = max_shift_c_min;
1699 }
1700
1701 // Now that barycenters and diagrams are updated in
1702 // PDBarycenter class, we import the results here.
1703 diagrams_c_min = barycenter_computer_min_[c].getCurrentBidders();
1704 centroids_with_price_min
1705 = barycenter_computer_min_[c].getCurrentBarycenter();
1706 int i = 0;
1707 for(int idx : clustering_[c]) {
1708 current_bidder_diagrams_min_[idx] = diagrams_c_min[i];
1709 centroids_with_price_min_[idx] = centroids_with_price_min[i];
1710 i++;
1711 }
1712
1713 GoodDiagram old_centroid = centroids_min_[c];
1714 centroids_min_[c] = centroidWithZeroPrices(
1715 centroids_with_price_min_[clustering_[c][0]]);
1716
1717 if(use_accelerated_) {
1718 wasserstein_shift
1719 += computeDistance(old_centroid, centroids_min_[c], 0.01);
1720 }
1721 }
1722
1723 if(do_sad_) {
1724 std::vector<std::vector<MatchingType>> all_matchings;
1725 total_cost = 0;
1726 std::vector<int> sizes;
1727
1728 std::vector<BidderDiagram> diagrams_c_min;
1729 if(barycenter_inputs_reset_flag) {
1730 for(int idx : clustering_[c]) {
1731 diagrams_c_min.emplace_back(current_bidder_diagrams_saddle_[idx]);
1732 }
1733 sizes.resize(diagrams_c_min.size());
1734 for(size_t i = 0; i < diagrams_c_min.size(); i++) {
1735 sizes[i] = diagrams_c_min[i].size();
1736 }
1737 barycenter_computer_sad_[c].setNumberOfInputs(diagrams_c_min.size());
1738 barycenter_computer_sad_[c].setCurrentBidders(diagrams_c_min);
1739 std::vector<GoodDiagram> barycenter_goods(clustering_[c].size());
1740 for(size_t i_diagram = 0; i_diagram < clustering_[c].size();
1741 i_diagram++) {
1742 barycenter_goods[i_diagram]
1743 = centroids_with_price_saddle_[clustering_[c][i_diagram]];
1744 }
1745 barycenter_computer_sad_[c].setCurrentBarycenter(barycenter_goods);
1746 all_matchings.resize(diagrams_c_min.size());
1747 } else {
1748 sizes.resize(barycenter_computer_sad_[c].getCurrentBidders().size());
1749 for(size_t i = 0;
1750 i < barycenter_computer_sad_[c].getCurrentBidders().size(); i++) {
1751 sizes[i]
1752 = barycenter_computer_sad_[c].getCurrentBidders().at(i).size();
1753 }
1754 all_matchings.resize(
1755 barycenter_computer_sad_[c].getCurrentBidders().size());
1756 diagrams_c_min.resize(
1757 barycenter_computer_sad_[c].getCurrentBidders().size());
1758 }
1759
1760 KDTreePair pair;
1761 bool use_kdt = false;
1762 if(!barycenter_computer_sad_[c].getCurrentBarycenter()[0].empty()) {
1763 pair = barycenter_computer_sad_[c].getKDTree();
1764 use_kdt = true;
1765 }
1766
1767 barycenter_computer_sad_[c].runMatching(
1768 &total_cost, epsilon_[1], sizes, *pair.first, pair.second,
1769 &(min_diag_price->at(1)), &(min_price->at(1)), &(all_matchings),
1770 use_kdt, only_matchings);
1771 for(size_t ii = 0; ii < all_matchings.size(); ii++) {
1772 all_matchings_per_type_and_cluster[c][1][ii].resize(
1773 all_matchings[ii].size());
1774 all_matchings_per_type_and_cluster[c][1][ii] = all_matchings[ii];
1775 }
1776 for(int ii = all_matchings.size(); ii < numberOfInputs_; ii++) {
1777 all_matchings_per_type_and_cluster[c][1][ii].resize(0);
1778 }
1779
1780 precision_sad
1781 = barycenter_computer_sad_[c].isPrecisionObjectiveMet(deltaLim_, 0);
1782 if(!only_matchings) {
1783 max_shift_c_sad
1784 = barycenter_computer_sad_[c].updateBarycenter(all_matchings);
1785 }
1786
1787 cost_sad_ += total_cost;
1788
1789 if(max_shift_c_sad > max_shift_vector[1]) {
1790 max_shift_vector[1] = max_shift_c_sad;
1791 }
1792
1793 // Now that barycenters and diagrams are updated in PDBarycenter class,
1794 // we import the results here.
1795 diagrams_c_min = barycenter_computer_sad_[c].getCurrentBidders();
1796 centroids_with_price_sad
1797 = barycenter_computer_sad_[c].getCurrentBarycenter();
1798 int i = 0;
1799 for(int idx : clustering_[c]) {
1800 current_bidder_diagrams_saddle_[idx] = diagrams_c_min[i];
1801 centroids_with_price_saddle_[idx] = centroids_with_price_sad[i];
1802 i++;
1803 }
1804 GoodDiagram old_centroid = centroids_saddle_[c];
1805 centroids_saddle_[c] = centroidWithZeroPrices(
1806 centroids_with_price_saddle_[clustering_[c][0]]);
1807 if(use_accelerated_)
1808 wasserstein_shift
1809 += computeDistance(old_centroid, centroids_saddle_[c], 0.01);
1810 }
1811
1812 if(do_max_) {
1813 std::vector<std::vector<MatchingType>> all_matchings;
1814 Timer time_preprocess_bary;
1815 total_cost = 0;
1816 std::vector<int> sizes;
1817 std::vector<BidderDiagram> diagrams_c_min;
1818 if(barycenter_inputs_reset_flag) {
1819 for(int idx : clustering_[c]) {
1820 diagrams_c_min.emplace_back(current_bidder_diagrams_max_[idx]);
1821 }
1822 sizes.resize(diagrams_c_min.size());
1823 for(size_t i = 0; i < diagrams_c_min.size(); i++) {
1824 sizes[i] = diagrams_c_min[i].size();
1825 }
1826 barycenter_computer_max_[c].setNumberOfInputs(diagrams_c_min.size());
1827 barycenter_computer_max_[c].setCurrentBidders(diagrams_c_min);
1828 std::vector<GoodDiagram> barycenter_goods(clustering_[c].size());
1829 for(size_t i_diagram = 0; i_diagram < clustering_[c].size();
1830 i_diagram++) {
1831 barycenter_goods[i_diagram]
1832 = centroids_with_price_max_[clustering_[c][i_diagram]];
1833 }
1834 barycenter_computer_max_[c].setCurrentBarycenter(barycenter_goods);
1835 all_matchings.resize(diagrams_c_min.size());
1836 } else {
1837 sizes.resize(barycenter_computer_max_[c].getCurrentBidders().size());
1838 for(size_t i = 0;
1839 i < barycenter_computer_max_[c].getCurrentBidders().size(); i++) {
1840 sizes[i]
1841 = barycenter_computer_max_[c].getCurrentBidders().at(i).size();
1842 }
1843
1844 diagrams_c_min.resize(
1845 barycenter_computer_max_[c].getCurrentBidders().size());
1846 all_matchings.resize(
1847 barycenter_computer_max_[c].getCurrentBidders().size());
1848 }
1849
1850 KDTreePair pair;
1851 bool use_kdt = false;
1852 if(!barycenter_computer_max_[c].getCurrentBarycenter()[0].empty()) {
1853 pair = barycenter_computer_max_[c].getKDTree();
1854 use_kdt = true;
1855 }
1856
1857 barycenter_computer_max_[c].runMatching(
1858 &total_cost, epsilon_[2], sizes, *pair.first, pair.second,
1859 &(min_diag_price->at(2)), &(min_price->at(2)), &(all_matchings),
1860 use_kdt, only_matchings);
1861 for(size_t ii = 0; ii < all_matchings.size(); ii++) {
1862 all_matchings_per_type_and_cluster[c][2][ii].resize(
1863 all_matchings[ii].size());
1864 all_matchings_per_type_and_cluster[c][2][ii] = all_matchings[ii];
1865 }
1866 for(int ii = all_matchings.size(); ii < numberOfInputs_; ii++) {
1867 all_matchings_per_type_and_cluster[c][2][ii].resize(0);
1868 }
1869 precision_max
1870 = barycenter_computer_max_[c].isPrecisionObjectiveMet(deltaLim_, 0);
1871
1872 cost_max_ += total_cost;
1873 Timer time_update;
1874 if(!only_matchings) {
1875 max_shift_c_max
1876 = barycenter_computer_max_[c].updateBarycenter(all_matchings);
1877 }
1878 if(max_shift_c_max > max_shift_vector[2]) {
1879 max_shift_vector[2] = max_shift_c_max;
1880 }
1881
1882 // Now that barycenters and diagrams are updated in PDBarycenter class,
1883 // we import the results here.
1884 diagrams_c_min = barycenter_computer_max_[c].getCurrentBidders();
1885 centroids_with_price_max
1886 = barycenter_computer_max_[c].getCurrentBarycenter();
1887 int i = 0;
1888 for(int idx : clustering_[c]) {
1889 current_bidder_diagrams_max_[idx] = diagrams_c_min[i];
1890 centroids_with_price_max_[idx] = centroids_with_price_max[i];
1891 i++;
1892 }
1893 GoodDiagram old_centroid = centroids_max_[c];
1894 centroids_max_[c] = centroidWithZeroPrices(
1895 centroids_with_price_max_[clustering_[c][0]]);
1896 if(use_accelerated_) {
1897 wasserstein_shift
1898 += computeDistance(old_centroid, centroids_max_[c], 0.01);
1899 }
1900 }
1901
1902 cost_ = cost_min_ + cost_sad_ + cost_max_;
1903 if(wasserstein_shift > max_wasserstein_shift) {
1904 max_wasserstein_shift = wasserstein_shift;
1905 }
1906 if(use_accelerated_) {
1907 for(int i = 0; i < numberOfInputs_; ++i) {
1908 // Step 5 of Accelerated KMeans: Update the lower bound on distance
1909 // thanks to the triangular inequality
1910 l_[i][c] = Geometry::pow(
1911 Geometry::pow(l_[i][c], 1. / wasserstein_)
1912 - Geometry::pow(wasserstein_shift, 1. / wasserstein_),
1913 wasserstein_);
1914 if(l_[i][c] < 0) {
1915 l_[i][c] = 0;
1916 }
1917 }
1918 for(int idx : clustering_[c]) {
1919 // Step 6, update the upper bound on the distance to the centroid
1920 // thanks to the triangle inequality
1921 u_[idx] = Geometry::pow(
1922 Geometry::pow(u_[idx], 1. / wasserstein_)
1923 + Geometry::pow(wasserstein_shift, 1. / wasserstein_),
1924 wasserstein_);
1925 r_[idx] = true;
1926 }
1927 }
1928 }
1929 }
1930 // Normally return max_shift, but it seems there is a bug
1931 // yielding max_shift > 100 * max_wasserstein_shift
1932 // which should logically not really happen...
1933 // This is supposed to be only a temporary patch...
1934 precision_min_ = precision_min;
1935 precision_sad_ = precision_sad;
1936 precision_max_ = precision_max;
1937 precision_criterion_ = precision_min && precision_sad && precision_max;
1938 barycenter_inputs_reset_flag = false;
1939 return max_shift_vector; // std::min(max_shift, max_wasserstein_shift);
1940}
1941
1943 for(int i = 0; i < numberOfInputs_; i++) {
1944 if(do_min_) {
1945 DiagramType *CTDiagram = &((*inputDiagramsMin_)[i]);
1946 BidderDiagram bidders;
1947 for(size_t j = 0; j < CTDiagram->size(); j++) {
1948 // Add bidder to bidders
1949 Bidder b((*CTDiagram)[j], j, lambda_);
1950
1951 b.setPositionInAuction(bidders.size());
1952 bidders.emplace_back(b);
1953 if(b.isDiagonal() || b.x_ == b.y_) {
1954 this->printMsg("Diagonal point in diagram", debug::Priority::DETAIL);
1955 }
1956 }
1957 bidder_diagrams_min_.emplace_back(bidders);
1958 current_bidder_diagrams_min_.emplace_back(BidderDiagram());
1959 centroids_with_price_min_.emplace_back(GoodDiagram());
1960 std::vector<int> ids(bidders.size());
1961 for(size_t j = 0; j < ids.size(); j++) {
1962 ids[j] = -1;
1963 }
1964 current_bidder_ids_min_.emplace_back(ids);
1965 }
1966
1967 if(do_sad_) {
1968 DiagramType *CTDiagram = &((*inputDiagramsSaddle_)[i]);
1969
1970 BidderDiagram bidders;
1971 for(size_t j = 0; j < CTDiagram->size(); j++) {
1972 // Add bidder to bidders
1973 Bidder b((*CTDiagram)[j], j, lambda_);
1974
1975 b.setPositionInAuction(bidders.size());
1976 bidders.emplace_back(b);
1977 if(b.isDiagonal() || b.x_ == b.y_) {
1978 this->printMsg("Diagonal point in diagram", debug::Priority::DETAIL);
1979 }
1980 }
1981 bidder_diagrams_saddle_.emplace_back(bidders);
1982 current_bidder_diagrams_saddle_.emplace_back(BidderDiagram());
1983 centroids_with_price_saddle_.emplace_back(GoodDiagram());
1984 std::vector<int> ids(bidders.size());
1985 for(size_t j = 0; j < ids.size(); j++) {
1986 ids[j] = -1;
1987 }
1988 current_bidder_ids_sad_.emplace_back(ids);
1989 }
1990
1991 if(do_max_) {
1992 DiagramType *CTDiagram = &((*inputDiagramsMax_)[i]);
1993
1994 BidderDiagram bidders;
1995 for(size_t j = 0; j < CTDiagram->size(); j++) {
1996 // Add bidder to bidders
1997 Bidder b((*CTDiagram)[j], j, lambda_);
1998
1999 b.setPositionInAuction(bidders.size());
2000 bidders.emplace_back(b);
2001 if(b.isDiagonal() || b.x_ == b.y_) {
2002 this->printMsg("Diagonal point in diagram", debug::Priority::DETAIL);
2003 }
2004 }
2005 bidder_diagrams_max_.emplace_back(bidders);
2006 current_bidder_diagrams_max_.emplace_back(BidderDiagram());
2007 centroids_with_price_max_.emplace_back(GoodDiagram());
2008 std::vector<int> ids(bidders.size());
2009 for(size_t j = 0; j < ids.size(); j++) {
2010 ids[j] = -1;
2011 }
2012 current_bidder_ids_max_.emplace_back(ids);
2013 }
2014 }
2015}
2016
2018 std::vector<double> &previous_min_persistence,
2019 std::vector<double> &min_persistence,
2020 std::vector<std::vector<double>> &initial_diagonal_prices,
2021 std::vector<std::vector<double>> &initial_off_diagonal_prices,
2022 std::vector<int> &min_points_to_add,
2023 bool add_points_to_barycenter,
2024 bool first_enrichment) {
2025
2026 std::vector<double> new_min_persistence = min_persistence;
2027
2028 if(!do_min_) {
2029 new_min_persistence[0] = previous_min_persistence[0];
2030 }
2031 if(!do_sad_) {
2032 new_min_persistence[1] = previous_min_persistence[1];
2033 }
2034 if(!do_max_) {
2035 new_min_persistence[2] = previous_min_persistence[2];
2036 }
2037
2038 // 1. Get size of the largest current diagram, deduce the maximal number of
2039 // points to append
2040 size_t max_diagram_size_min{};
2041 size_t max_diagram_size_sad{};
2042 size_t max_diagram_size_max{};
2043 if(do_min_) {
2044 for(int i = 0; i < numberOfInputs_; i++) {
2045 if(current_bidder_diagrams_min_[i].size() > max_diagram_size_min) {
2046 max_diagram_size_min = current_bidder_diagrams_min_[i].size();
2047 }
2048 }
2049 }
2050 if(do_sad_) {
2051 for(int i = 0; i < numberOfInputs_; i++) {
2052 if(current_bidder_diagrams_saddle_[i].size() > max_diagram_size_sad) {
2053 max_diagram_size_sad = current_bidder_diagrams_saddle_[i].size();
2054 }
2055 }
2056 }
2057 if(do_max_) {
2058 for(int i = 0; i < numberOfInputs_; i++) {
2059 if(current_bidder_diagrams_max_[i].size() > max_diagram_size_max) {
2060 max_diagram_size_max = current_bidder_diagrams_max_[i].size();
2061 }
2062 }
2063 }
2064 int max_points_to_add_min
2065 = std::max(min_points_to_add[0],
2066 min_points_to_add[0] + (int)(max_diagram_size_min / 10));
2067 int max_points_to_add_sad
2068 = std::max(min_points_to_add[1],
2069 min_points_to_add[1] + (int)(max_diagram_size_sad / 10));
2070 int max_points_to_add_max
2071 = std::max(min_points_to_add[2],
2072 min_points_to_add[2] + (int)(max_diagram_size_max / 10));
2073
2074 // 2. Get which points can be added, deduce the new minimal persistence
2075 std::vector<std::vector<int>> candidates_to_be_added_min(numberOfInputs_);
2076 std::vector<std::vector<int>> candidates_to_be_added_sad(numberOfInputs_);
2077 std::vector<std::vector<int>> candidates_to_be_added_max(numberOfInputs_);
2078 std::vector<std::vector<int>> idx_min(numberOfInputs_);
2079 std::vector<std::vector<int>> idx_sad(numberOfInputs_);
2080 std::vector<std::vector<int>> idx_max(numberOfInputs_);
2081
2082 if(do_min_) {
2083 for(int i = 0; i < numberOfInputs_; i++) {
2084 std::vector<double> persistences;
2085 for(size_t j = 0; j < bidder_diagrams_min_[i].size(); j++) {
2086 Bidder b = bidder_diagrams_min_[i].at(j);
2087 double persistence = b.getPersistence();
2088 if(persistence >= min_persistence[0]
2089 && persistence <= previous_min_persistence[0]) {
2090 candidates_to_be_added_min[i].emplace_back(j);
2091 idx_min[i].emplace_back(idx_min[i].size());
2092 persistences.emplace_back(persistence);
2093 }
2094 }
2095 sort(
2096 idx_min[i].begin(), idx_min[i].end(), [&persistences](int &a, int &b) {
2097 return ((persistences[a] > persistences[b])
2098 || ((persistences[a] == persistences[b]) && (a > b)));
2099 });
2100 int size = candidates_to_be_added_min[i].size();
2101 if(size >= max_points_to_add_min) {
2102 double last_persistence_added_min
2103 = persistences[idx_min[i][max_points_to_add_min - 1]];
2104 if(first_enrichment) { // a minima min_point_to_add (=max_point_to_add)
2105 // added per diagram
2106 if(i == 0) {
2107 new_min_persistence[0] = last_persistence_added_min;
2108 } else {
2109 if(last_persistence_added_min < new_min_persistence[0])
2110 new_min_persistence[0] = last_persistence_added_min;
2111 }
2112 } else { // a maxima max_point_to_add added per diagram
2113 if(last_persistence_added_min > new_min_persistence[0]) {
2114 new_min_persistence[0] = last_persistence_added_min;
2115 }
2116 }
2117 }
2118 }
2119 }
2120
2121 if(do_sad_) {
2122 for(int i = 0; i < numberOfInputs_; i++) {
2123 std::vector<double> persistences;
2124 for(size_t j = 0; j < bidder_diagrams_saddle_[i].size(); j++) {
2125 Bidder b = bidder_diagrams_saddle_[i].at(j);
2126 double persistence = b.getPersistence();
2127 if(persistence >= min_persistence[1]
2128 && persistence <= previous_min_persistence[1]) {
2129 candidates_to_be_added_sad[i].emplace_back(j);
2130 idx_sad[i].emplace_back(idx_sad[i].size());
2131 persistences.emplace_back(persistence);
2132 }
2133 }
2134 sort(
2135 idx_sad[i].begin(), idx_sad[i].end(), [&persistences](int &a, int &b) {
2136 return ((persistences[a] > persistences[b])
2137 || ((persistences[a] == persistences[b]) && (a > b)));
2138 });
2139 int size = candidates_to_be_added_sad[i].size();
2140 if(size >= max_points_to_add_sad) {
2141 double last_persistence_added_sad
2142 = persistences[idx_sad[i][max_points_to_add_sad - 1]];
2143 if(first_enrichment) { // a minima min_point_to_add (=max_point_to_add)
2144 // added per diagram
2145 if(i == 0) {
2146 new_min_persistence[1] = last_persistence_added_sad;
2147 } else {
2148 if(last_persistence_added_sad < new_min_persistence[1])
2149 new_min_persistence[1] = last_persistence_added_sad;
2150 }
2151 } else { // a maxima max_point_to_add added per diagram
2152 if(last_persistence_added_sad > new_min_persistence[1]) {
2153 new_min_persistence[1] = last_persistence_added_sad;
2154 }
2155 }
2156 }
2157 }
2158 }
2159 if(do_max_) {
2160 for(int i = 0; i < numberOfInputs_; i++) {
2161 std::vector<double> persistences;
2162 for(size_t j = 0; j < bidder_diagrams_max_[i].size(); j++) {
2163 Bidder b = bidder_diagrams_max_[i].at(j);
2164 double persistence = b.getPersistence();
2165 if(persistence >= min_persistence[2]
2166 && persistence <= previous_min_persistence[2]) {
2167 candidates_to_be_added_max[i].emplace_back(j);
2168 idx_max[i].emplace_back(idx_max[i].size());
2169 persistences.emplace_back(persistence);
2170 }
2171 }
2172 sort(
2173 idx_max[i].begin(), idx_max[i].end(), [&persistences](int &a, int &b) {
2174 return ((persistences[a] > persistences[b])
2175 || ((persistences[a] == persistences[b]) && (a > b)));
2176 });
2177 int size = candidates_to_be_added_max[i].size();
2178 if(size >= max_points_to_add_max) {
2179 double last_persistence_added_max
2180 = persistences[idx_max[i][max_points_to_add_max - 1]];
2181 if(first_enrichment) { // a minima min_point_to_add (=max_point_to_add)
2182 // added per diagram
2183 if(i == 0) {
2184 new_min_persistence[2] = last_persistence_added_max;
2185 } else {
2186 if(last_persistence_added_max < new_min_persistence[2])
2187 new_min_persistence[2] = last_persistence_added_max;
2188 }
2189 } else { // a maxima max_point_to_add added per diagram
2190 if(last_persistence_added_max > new_min_persistence[2]) {
2191 new_min_persistence[2] = last_persistence_added_max;
2192 }
2193 }
2194 }
2195 }
2196 }
2197
2198 // 3. Add the points to the current diagrams
2199 if(do_min_) {
2200 int compteur_for_adding_points = 0;
2201 for(int i = 0; i < numberOfInputs_; i++) {
2202 int size = candidates_to_be_added_min[i].size();
2203 for(int j = 0; j < std::min(max_points_to_add_min, size); j++) {
2204 Bidder b = bidder_diagrams_min_[i].at(
2205 candidates_to_be_added_min[i][idx_min[i][j]]);
2206 double persistence = b.getPersistence();
2207 if(persistence >= new_min_persistence[0]) {
2208 b.id_ = current_bidder_diagrams_min_[i].size();
2209 b.setPositionInAuction(current_bidder_diagrams_min_[i].size());
2210 b.setDiagonalPrice(initial_diagonal_prices[0][i]);
2211 current_bidder_diagrams_min_[i].emplace_back(b);
2212 current_bidder_ids_min_[i]
2213 [candidates_to_be_added_min[i][idx_min[i][j]]]
2214 = current_bidder_diagrams_min_[i].size() - 1;
2215
2216 if(use_accelerated_ && n_iterations_ > 0) {
2217 for(int c = 0; c < k_; ++c) {
2218 // Step 5 of Accelerated KMeans: Update the lower bound on
2219 // distance thanks to the triangular inequality
2220 l_[i][c]
2221 = Geometry::pow(Geometry::pow(l_[i][c], 1. / wasserstein_)
2222 - persistence / sqrt(2),
2223 wasserstein_);
2224 if(l_[i][c] < 0) {
2225 l_[i][c] = 0;
2226 }
2227 }
2228 // Step 6, update the upper bound on the distance to the centroid
2229 // thanks to the triangle inequality
2230 u_[i] = Geometry::pow(
2231 Geometry::pow(u_[i], 1. / wasserstein_) + persistence / sqrt(2),
2232 wasserstein_);
2233 r_[i] = true;
2234 }
2235 int to_be_added_to_barycenter
2236 = deterministic_ ? compteur_for_adding_points % numberOfInputs_
2237 : rand() % numberOfInputs_;
2238 if(to_be_added_to_barycenter == 0 && add_points_to_barycenter) {
2239 for(int k = 0; k < numberOfInputs_; k++) {
2240 if(inv_clustering_[i] == inv_clustering_[k]) {
2241 Good g = Good(
2242 b.x_, b.y_, false, centroids_with_price_min_[k].size());
2243 g.setPrice(initial_off_diagonal_prices[0][k]);
2245 b.coords_[0], b.coords_[1], b.coords_[2]);
2246 centroids_with_price_min_[k].emplace_back(g);
2247 }
2248 }
2249 Good g = Good(
2250 b.x_, b.y_, false, centroids_min_[inv_clustering_[i]].size());
2251 g.SetCriticalCoordinates(b.coords_[0], b.coords_[1], b.coords_[2]);
2252 centroids_min_[inv_clustering_[i]].emplace_back(g);
2253 }
2254 }
2255 compteur_for_adding_points++;
2256 }
2257 }
2258 }
2259 if(do_sad_) {
2260 int compteur_for_adding_points = 0;
2261 for(int i = 0; i < numberOfInputs_; i++) {
2262 int size = candidates_to_be_added_sad[i].size();
2263 for(int j = 0; j < std::min(max_points_to_add_sad, size); j++) {
2264 Bidder b = bidder_diagrams_saddle_[i].at(
2265 candidates_to_be_added_sad[i][idx_sad[i][j]]);
2266 double persistence = b.getPersistence();
2267 if(persistence >= new_min_persistence[1]) {
2268 b.id_ = current_bidder_diagrams_saddle_[i].size();
2269 b.setPositionInAuction(current_bidder_diagrams_saddle_[i].size());
2270 b.setDiagonalPrice(initial_diagonal_prices[1][i]);
2271 current_bidder_diagrams_saddle_[i].emplace_back(b);
2272 current_bidder_ids_sad_[i]
2273 [candidates_to_be_added_sad[i][idx_sad[i][j]]]
2274 = current_bidder_diagrams_saddle_[i].size() - 1;
2275
2276 if(use_accelerated_ && n_iterations_ > 0) {
2277 for(int c = 0; c < k_; ++c) {
2278 // Step 5 of Accelerated KMeans: Update the lower bound on
2279 // distance thanks to the triangular inequality
2280 l_[i][c]
2281 = Geometry::pow(Geometry::pow(l_[i][c], 1. / wasserstein_)
2282 - persistence / sqrt(2),
2283 wasserstein_);
2284 if(l_[i][c] < 0) {
2285 l_[i][c] = 0;
2286 }
2287 }
2288 // Step 6, update the upper bound on the distance to the centroid
2289 // thanks to the triangle inequality
2290 u_[i] = Geometry::pow(
2291 Geometry::pow(u_[i], 1. / wasserstein_) + persistence / sqrt(2),
2292 wasserstein_);
2293 r_[i] = true;
2294 }
2295 int to_be_added_to_barycenter
2296 = deterministic_ ? compteur_for_adding_points % numberOfInputs_
2297 : rand() % numberOfInputs_;
2298 if(to_be_added_to_barycenter == 0 && add_points_to_barycenter) {
2299 for(int k = 0; k < numberOfInputs_; k++) {
2300 if(inv_clustering_[i] == inv_clustering_[k]) {
2301 Good g = Good(
2302 b.x_, b.y_, false, centroids_with_price_saddle_[k].size());
2303 g.setPrice(initial_off_diagonal_prices[1][k]);
2305 b.coords_[0], b.coords_[1], b.coords_[2]);
2306 centroids_with_price_saddle_[k].emplace_back(g);
2307 }
2308 }
2309 }
2310 }
2311 compteur_for_adding_points++;
2312 }
2313 }
2314 }
2315 if(do_max_) {
2316 int compteur_for_adding_points = 0;
2317 for(int i = 0; i < numberOfInputs_; i++) {
2318 int size = candidates_to_be_added_max[i].size();
2319 for(int j = 0; j < std::min(max_points_to_add_max, size); j++) {
2320 Bidder b = bidder_diagrams_max_[i].at(
2321 candidates_to_be_added_max[i][idx_max[i][j]]);
2322 double persistence = b.getPersistence();
2323 if(persistence >= new_min_persistence[2]) {
2324 b.id_ = current_bidder_diagrams_max_[i].size();
2325 b.setPositionInAuction(current_bidder_diagrams_max_[i].size());
2326 b.setDiagonalPrice(initial_diagonal_prices[2][i]);
2327 current_bidder_diagrams_max_[i].emplace_back(b);
2328 current_bidder_ids_max_[i]
2329 [candidates_to_be_added_max[i][idx_max[i][j]]]
2330 = current_bidder_diagrams_max_[i].size() - 1;
2331
2332 if(use_accelerated_ && n_iterations_ > 0) {
2333 for(int c = 0; c < k_; ++c) {
2334 // Step 5 of Accelerated KMeans: Update the lower bound on
2335 // distance thanks to the triangular inequality
2336 l_[i][c]
2337 = Geometry::pow(Geometry::pow(l_[i][c], 1. / wasserstein_)
2338 - persistence / sqrt(2),
2339 wasserstein_);
2340 if(l_[i][c] < 0) {
2341 l_[i][c] = 0;
2342 }
2343 }
2344 // Step 6, update the upper bound on the distance to the centroid
2345 // thanks to the triangle inequality
2346 u_[i] = Geometry::pow(
2347 Geometry::pow(u_[i], 1. / wasserstein_) + persistence / sqrt(2),
2348 wasserstein_);
2349 r_[i] = true;
2350 }
2351 int to_be_added_to_barycenter
2352 = deterministic_ ? compteur_for_adding_points % numberOfInputs_
2353 : rand() % numberOfInputs_;
2354 if(to_be_added_to_barycenter == 0 && add_points_to_barycenter) {
2355 for(int k = 0; k < numberOfInputs_; k++) {
2356 if(inv_clustering_[i] == inv_clustering_[k]) {
2357 Good g = Good(
2358 b.x_, b.y_, false, centroids_with_price_max_[k].size());
2359 g.setPrice(initial_off_diagonal_prices[2][k]);
2361 b.coords_[0], b.coords_[1], b.coords_[2]);
2362 centroids_with_price_max_[k].emplace_back(g);
2363 }
2364 }
2365 Good g = Good(
2366 b.x_, b.y_, false, centroids_max_[inv_clustering_[i]].size());
2367 g.SetCriticalCoordinates(b.coords_[0], b.coords_[1], b.coords_[2]);
2368 centroids_max_[inv_clustering_[i]].emplace_back(g);
2369 }
2370 }
2371 compteur_for_adding_points++;
2372 }
2373 }
2374 }
2375
2376 return new_min_persistence;
2377}
2378
2380 std::vector<double> & /*min_persistence*/) {
2381
2382 if(do_min_) {
2383 barycenter_computer_min_.resize(k_);
2384 for(int c = 0; c < k_; c++) {
2385 std::vector<BidderDiagram> diagrams_c;
2386 for(int idx : clustering_[c]) {
2387 diagrams_c.emplace_back(current_bidder_diagrams_min_[idx]);
2388 }
2389 barycenter_computer_min_[c] = {};
2390 barycenter_computer_min_[c].setThreadNumber(threadNumber_);
2391 barycenter_computer_min_[c].setWasserstein(wasserstein_);
2392 barycenter_computer_min_[c].setDiagramType(0);
2393 barycenter_computer_min_[c].setUseProgressive(false);
2394 barycenter_computer_min_[c].setDeterministic(true);
2395 barycenter_computer_min_[c].setGeometricalFactor(geometrical_factor_);
2396 barycenter_computer_min_[c].setDebugLevel(debugLevel_);
2397 barycenter_computer_min_[c].setNumberOfInputs(diagrams_c.size());
2398 barycenter_computer_min_[c].setCurrentBidders(diagrams_c);
2399 }
2400 }
2401 if(do_sad_) {
2402 barycenter_computer_sad_.resize(k_);
2403 for(int c = 0; c < k_; c++) {
2404 std::vector<BidderDiagram> diagrams_c;
2405 for(int idx : clustering_[c]) {
2406 diagrams_c.emplace_back(current_bidder_diagrams_saddle_[idx]);
2407 }
2408 barycenter_computer_sad_[c] = {};
2409 barycenter_computer_sad_[c].setThreadNumber(threadNumber_);
2410 barycenter_computer_sad_[c].setWasserstein(wasserstein_);
2411 barycenter_computer_sad_[c].setDiagramType(1);
2412 barycenter_computer_sad_[c].setUseProgressive(false);
2413 barycenter_computer_sad_[c].setDeterministic(true);
2414 barycenter_computer_sad_[c].setGeometricalFactor(geometrical_factor_);
2415 barycenter_computer_sad_[c].setDebugLevel(debugLevel_);
2416 barycenter_computer_sad_[c].setNumberOfInputs(diagrams_c.size());
2417 barycenter_computer_sad_[c].setCurrentBidders(diagrams_c);
2418
2419 std::vector<GoodDiagram> barycenter_goods(clustering_[c].size());
2420 for(size_t i_diagram = 0; i_diagram < clustering_[c].size();
2421 i_diagram++) {
2422 barycenter_goods[i_diagram]
2423 = centroids_with_price_saddle_[clustering_[c][i_diagram]];
2424 }
2425 barycenter_computer_sad_[c].setCurrentBarycenter(barycenter_goods);
2426 }
2427 }
2428 if(do_max_) {
2429 barycenter_computer_max_.resize(k_);
2430 for(int c = 0; c < k_; c++) {
2431 std::vector<BidderDiagram> diagrams_c;
2432 for(int idx : clustering_[c]) {
2433 diagrams_c.emplace_back(current_bidder_diagrams_max_[idx]);
2434 }
2435 barycenter_computer_max_[c] = {};
2436 barycenter_computer_max_[c].setDiagrams(inputDiagramsMax_);
2437 barycenter_computer_max_[c].setThreadNumber(threadNumber_);
2438 barycenter_computer_max_[c].setWasserstein(wasserstein_);
2439 barycenter_computer_max_[c].setDiagramType(2);
2440 barycenter_computer_max_[c].setUseProgressive(false);
2441 barycenter_computer_max_[c].setDeterministic(true);
2442 barycenter_computer_max_[c].setGeometricalFactor(geometrical_factor_);
2443 barycenter_computer_max_[c].setDebugLevel(debugLevel_);
2444 barycenter_computer_max_[c].setNumberOfInputs(diagrams_c.size());
2445 barycenter_computer_max_[c].setCurrentBidders(diagrams_c);
2446
2447 std::vector<GoodDiagram> barycenter_goods(clustering_[c].size());
2448 for(size_t i_diagram = 0; i_diagram < clustering_[c].size();
2449 i_diagram++) {
2450 barycenter_goods[i_diagram]
2451 = centroids_with_price_max_[clustering_[c][i_diagram]];
2452 }
2453 barycenter_computer_max_[c].setCurrentBarycenter(barycenter_goods);
2454 }
2455 }
2456}
2457
2459 std::ofstream ufile("u_vec.txt");
2460 std::ofstream lfile("l_mat.txt");
2461 std::ofstream approx_file("a_mat.txt");
2462 if(ufile.is_open() && lfile.is_open()) {
2463 for(int i = 0; i < numberOfInputs_; i++) {
2464 ufile << u_[i] << " ";
2465 for(int j = 0; j < k_; j++) {
2466 lfile << l_[i][j] << " ";
2467 }
2468 lfile << "\n";
2469 }
2470 }
2471
2472 for(int c = 0; c < k_; c++) {
2473 for(int i : clustering_[c]) {
2474 approx_file << (u_[i] + l_[i][c]) / 2 << " ";
2475 }
2476 approx_file << "\n";
2477 }
2478 lfile.close();
2479 ufile.close();
2480 approx_file.close();
2481}
2482
2484 std::cout << "Computing real distances to every clusters" << std::endl;
2485 std::ofstream file("a_real_mat.txt");
2486 if(file.is_open()) {
2487 for(int c = 0; c < k_; c++) {
2488 for(int i : clustering_[c]) {
2489 file << distanceToCentroid_[i] << " ";
2490 }
2491 file << "\n";
2492 }
2493 file.close();
2494 } else {
2495 std::cout << "file not open" << std::endl;
2496 }
2497}
2498
2500 std::ofstream file(
2501 "prices_evolution.txt", std::ofstream::out | std::ofstream::app);
2502 if(file.is_open()) {
2503 file << "\nITERATION " << iteration << "\n" << std::endl;
2504 for(int i = 0; i < k_; i++) {
2505 file << "\ncentroid " << i << std::endl;
2506
2507 for(size_t j = 0; j < centroids_with_price_max_[i].size(); j++) {
2508 Good g = centroids_with_price_max_[i].at(j);
2509 file << g.getPrice() << " ";
2510 }
2511 }
2512 }
2513 file.close();
2514}
2515
2517 double total_real_cost_min = 0;
2518 double total_real_cost_max = 0;
2519 double total_real_cost_sad = 0;
2520 double sq_distance;
2521 if(original_dos[0]) {
2522 for(int c = 0; c < k_; c++) {
2523 double real_cost_cluster = 0;
2524 for(int i = 0; i < numberOfInputs_; i++) {
2525 GoodDiagram current_barycenter
2526 = centroidWithZeroPrices(centroids_min_[c]);
2527 sq_distance
2528 = computeDistance(bidder_diagrams_min_[i], current_barycenter, 0.01);
2529 real_cost_cluster += sq_distance;
2530 }
2531 total_real_cost_min += real_cost_cluster;
2532 }
2533 }
2534 if(original_dos[1]) {
2535 for(int c = 0; c < k_; c++) {
2536 double real_cost_cluster = 0;
2537 for(int i = 0; i < numberOfInputs_; i++) {
2538 GoodDiagram current_barycenter
2539 = centroidWithZeroPrices(centroids_saddle_[c]);
2540 sq_distance = computeDistance(
2541 bidder_diagrams_saddle_[i], current_barycenter, 0.01);
2542 real_cost_cluster += sq_distance;
2543 }
2544 total_real_cost_sad += real_cost_cluster;
2545 }
2546 }
2547 if(original_dos[2]) {
2548 for(int c = 0; c < k_; c++) {
2549 double real_cost_cluster = 0;
2550 for(int i = 0; i < numberOfInputs_; i++) {
2551 GoodDiagram current_barycenter
2552 = centroidWithZeroPrices(centroids_max_[c]);
2553 sq_distance
2554 = computeDistance(bidder_diagrams_max_[i], current_barycenter, 0.01);
2555 real_cost_cluster += sq_distance;
2556 }
2557 total_real_cost_max += real_cost_cluster;
2558 }
2559 }
2560 return total_real_cost_min + total_real_cost_sad + total_real_cost_max;
2561}
2562
2564 std::vector<std::vector<std::vector<std::vector<MatchingType>>>>
2565 &all_matchings_per_type_and_cluster) {
2566
2567 if(do_min_) {
2568 computeBarycenterForTwo(
2569 all_matchings_per_type_and_cluster[0][0], current_bidder_ids_min_,
2570 current_bidder_diagrams_min_, bidder_diagrams_min_, centroids_min_[0]);
2571 }
2572 if(do_sad_) {
2573 computeBarycenterForTwo(all_matchings_per_type_and_cluster[0][1],
2574 current_bidder_ids_sad_,
2575 current_bidder_diagrams_saddle_,
2576 bidder_diagrams_saddle_, centroids_saddle_[0]);
2577 }
2578 if(do_max_) {
2579 computeBarycenterForTwo(
2580 all_matchings_per_type_and_cluster[0][2], current_bidder_ids_max_,
2581 current_bidder_diagrams_max_, bidder_diagrams_max_, centroids_max_[0]);
2582 }
2583}
2584
2586 std::vector<std::vector<MatchingType>> &matchings,
2587 std::vector<std::vector<int>> &bidders_ids,
2588 std::vector<BidderDiagram> &current_bidder_diagrams,
2589 std::vector<BidderDiagram> &bidder_diagrams,
2590 GoodDiagram &barycenter) {
2591
2592 auto &matchings0 = matchings[0];
2593 auto &diagram0 = bidder_diagrams[0];
2594 auto &current_diagram0 = current_bidder_diagrams[0];
2595 auto &ids0 = bidders_ids[0];
2596
2597 auto &matchings1 = matchings[1];
2598 auto &diagram1 = bidder_diagrams[1];
2599 auto &current_diagram1 = current_bidder_diagrams[1];
2600 auto &ids1 = bidders_ids[1];
2601
2602 std::vector<int> new_to_old_id(current_diagram1.size());
2603 // 1. Invert the current_bidder_ids_ vector
2604 for(size_t j = 0; j < ids1.size(); j++) {
2605 int new_id = ids1[j];
2606 if(new_id >= 0) {
2607 new_to_old_id[new_id] = j;
2608 }
2609 }
2610
2611 std::vector<MatchingType> matching_to_add(0);
2612 std::vector<MatchingType> matching_to_add2(0);
2613
2614 for(size_t i = 0; i < matchings1.size(); i++) {
2615 MatchingType &t = matchings1[i];
2616 int bidderId = std::get<0>(t);
2617 int goodId = std::get<1>(t);
2618
2619 if(bidderId >= 0) {
2620 Bidder b = diagram1.at(new_to_old_id[bidderId]);
2621 double bx = b.x_;
2622 double by = b.y_;
2623 if(goodId >= 0) {
2624 Good g = barycenter.at(goodId);
2625 double gx = g.x_;
2626 double gy = g.y_;
2627 barycenter.at(goodId).x_ = (bx + gx) / 2;
2628 barycenter.at(goodId).y_ = (by + gy) / 2;
2629 // divide by 4 in order to display the cost of the half matching
2630 // i.e. the cost of matching to the barycenter
2631 std::get<2>(t) /= 4;
2632
2633 } else {
2634 double gx = (bx + by) / 2;
2635 double gy = (bx + by) / 2;
2636 gx = (gx + bx) / 2;
2637 gy = (gy + by) / 2;
2638 double cost = Geometry::pow((gx - bx), wasserstein_)
2639 + Geometry::pow((gy - by), wasserstein_);
2640 MatchingType t2 = std::make_tuple(bidderId, barycenter.size(), cost);
2641 MatchingType t3 = std::make_tuple(-1, barycenter.size(), cost);
2642 Good g = Good(gx, gy, false, barycenter.size());
2643 // g.SetCriticalCoordinates(b.coords_[0], b.coords_[1], b.coords_[2]);
2644 barycenter.emplace_back(g);
2645 matching_to_add.emplace_back(t2);
2646 matching_to_add2.emplace_back(t3);
2647 }
2648 } else {
2649 if(goodId >= 0) {
2650 Good g = barycenter.at(goodId);
2651 double gx = (g.x_ + g.y_) / 2;
2652 double gy = (g.x_ + g.y_) / 2;
2653 barycenter.at(goodId).x_ = (gx + g.x_) / 2;
2654 barycenter.at(goodId).y_ = (gy + g.y_) / 2;
2655 std::get<2>(t) /= 4;
2656 }
2657 }
2658 }
2659 for(size_t j = 0; j < matching_to_add.size(); j++) {
2660 matchings1.emplace_back(matching_to_add[j]);
2661 }
2662 for(size_t j = 0; j < matching_to_add2.size(); j++) {
2663 matchings0.emplace_back(matching_to_add2[j]);
2664 }
2665
2666 // correct the costs of matchings in diagram 0
2667 // costs are initially allzeros because the barycenter is identical
2668 // to diagram 0
2669 std::vector<int> new_to_old_id2(current_diagram0.size());
2670 // 1. Invert the current_bidder_ids_ vector
2671 for(size_t j = 0; j < ids0.size(); j++) {
2672 int new_id = ids0[j];
2673 if(new_id >= 0) {
2674 new_to_old_id2[new_id] = j;
2675 }
2676 }
2677
2678 for(size_t i = 0; i < matchings0.size(); i++) {
2679 MatchingType &t = matchings0[i];
2680 int bidderId = std::get<0>(t);
2681 int goodId = std::get<1>(t);
2682
2683 if(bidderId >= 0 and goodId >= 0) {
2684 Bidder b = diagram0.at(new_to_old_id2[bidderId]);
2685 double bx = b.x_;
2686 double by = b.y_;
2687 Good g = barycenter.at(goodId);
2688 double gx = g.x_;
2689 double gy = g.y_;
2690 double cost = Geometry::pow((gx - bx), wasserstein_)
2691 + Geometry::pow((gy - by), wasserstein_);
2692 std::get<2>(t) = cost;
2693 }
2694 }
2695}
void setDiagonalPrice(const double price)
void setPositionInAuction(const int pos)
int debugLevel_
Definition: Debug.h:379
int printWrn(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition: Debug.h:159
int printMsg(const std::string &msg, const debug::Priority &priority=debug::Priority::INFO, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cout) const
Definition: Debug.h:118
void setPrice(const double price)
std::pair< typename KDT::KDTreeRoot, typename KDT::KDTreeMap > KDTreePair
Definition: PDBarycenter.h:58
std::vector< std::vector< int > > old_clustering_
Definition: PDClustering.h:308
void initializeAcceleratedKMeans()
std::vector< double > enrichCurrentBidderDiagrams(std::vector< double > &previous_min_persistence, std::vector< double > &min_persistence, std::vector< std::vector< double > > &initial_diagonal_prices, std::vector< std::vector< double > > &initial_off_diagonal_prices, std::vector< int > &min_points_to_add, bool add_points_to_barycenter, bool first_enrichment)
std::vector< int > inv_clustering_
Definition: PDClustering.h:309
void computeDistanceToCentroid()
double getLessPersistent(int type=-1)
void computeBarycenterForTwo(std::vector< std::vector< MatchingType > > &matchings, std::vector< std::vector< int > > &bidders_ids, std::vector< BidderDiagram > &current_bidder_diagrams, std::vector< BidderDiagram > &bidder_diagrams, GoodDiagram &barycenter)
BidderDiagram centroidToDiagram(const GoodDiagram &centroid)
void initializeCentroidsKMeanspp()
double computeDistance(const BidderDiagram &D1, const BidderDiagram &D2, const double delta_lim)
void printRealDistancesToFile()
bool barycenter_inputs_reset_flag
Definition: PDClustering.h:246
void getCentroidDistanceMatrix()
std::vector< double > updateCentroidsPosition(std::vector< std::vector< double > > *min_price, std::vector< std::vector< double > > *min_diag_price, std::vector< std::vector< std::vector< std::vector< MatchingType > > > > &all_matchings, int only_matchings)
void printPricesToFile(int)
BidderDiagram diagramWithZeroPrices(const BidderDiagram &diagram)
double getMostPersistent(int type=-1)
std::vector< std::vector< int > > clustering_
Definition: PDClustering.h:307
std::vector< GoodDiagram > centroids_min_
Definition: PDClustering.h:292
GoodDiagram diagramToCentroid(const BidderDiagram &diagram)
void acceleratedUpdateClusters()
std::array< double, 3 > epsilon_
Definition: PDClustering.h:274
std::vector< GoodDiagram > centroids_saddle_
Definition: PDClustering.h:298
std::vector< std::vector< double > > getMinDiagonalPrices()
void computeBarycenterForTwoGlobal(std::vector< std::vector< std::vector< std::vector< MatchingType > > > > &)
std::vector< std::vector< int > > centroids_sizes_
Definition: PDClustering.h:311
void initializeEmptyClusters()
std::vector< std::vector< double > > getMinPrices()
void printMatchings(std::vector< std::vector< std::vector< MatchingType > > >)
std::vector< std::vector< double > > getDistanceMatrix()
void initializeBarycenterComputers(std::vector< double > &min_persistence)
GoodDiagram centroidWithZeroPrices(const GoodDiagram &centroid)
void resetDosToOriginalValues()
Definition: PDClustering.h:119
std::vector< GoodDiagram > centroids_max_
Definition: PDClustering.h:304
std::vector< int > execute(std::vector< DiagramType > &final_centroids, std::vector< std::vector< std::vector< std::vector< MatchingType > > > > &all_matchings)
std::array< bool, 3 > original_dos
Definition: PDClustering.h:287
std::vector< std::vector< int > > get_centroids_sizes()
void correctMatchings(std::vector< std::vector< std::vector< std::vector< MatchingType > > > > &previous_matchings)
void SetCriticalCoordinates(const float coords_x, const float coords_y, const float coords_z)
std::array< float, 3 > GetCriticalCoordinates() const
double run(std::vector< MatchingType > &matchings, const int kdt_index=0)
void BuildAuctionDiagrams(const BidderDiagram &BD, const GoodDiagram &GD)
double getElapsedTime()
Definition: Timer.h:15
T1 pow(const T1 val, const T2 n)
Definition: Geometry.h:411
std::vector< Bidder > BidderDiagram
CriticalType
default value for critical index
Definition: DataTypes.h:80
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)