TTK
Loading...
Searching...
No Matches
BottleneckDistance.cpp
Go to the documentation of this file.
1#include <AssignmentMunkres.h>
3#include <GabowTarjan.h>
4#include <Geometry.h>
5
7 this->setDebugMsgPrefix("BottleneckDistance");
8}
9
11 const ttk::DiagramType &diag1,
12 std::vector<MatchingType> &matchings) {
13 Timer t;
14
15 if(diag0.empty() || diag1.empty()) {
16 this->printErr("Empty input diagrams");
17 return -1;
18 }
19
20 bool fromParaView = this->PVAlgorithm >= 0;
21 if(fromParaView) {
22 switch(this->PVAlgorithm) {
23 case 0:
24 this->printMsg("Solving with the TTK approach");
25 this->computeBottleneck(diag0, diag1, matchings);
26 break;
27 case 1: {
28 this->printMsg("Solving with the legacy Dionysus exact approach.");
29 this->printErr("Not supported");
30 } break;
31 case 2: {
32 this->printMsg(
33 "Solving with the approximate Dionysus geometric approach.");
34 this->printErr("Not supported");
35 } break;
36 case 3: {
37 this->printMsg("Solving with the parallel TTK approach");
38 this->printErr("Not supported");
39 } break;
40 case 4: {
41 this->printMsg("Benchmarking");
42 this->printErr("Not supported");
43 } break;
44 default: {
45 this->printErr("You must specify a valid assignment algorithm.");
46 }
47 }
48
49 } else {
50 switch(str2int(this->DistanceAlgorithm.c_str())) {
51 case str2int("0"):
52 case str2int("ttk"):
53 this->printMsg("Solving with the TTK approach");
54 this->computeBottleneck(diag0, diag1, matchings);
55 break;
56 case str2int("1"):
57 case str2int("legacy"): {
58 this->printMsg("Solving with the legacy Dionysus exact approach.");
59 this->printErr("Not supported");
60 } break;
61 case str2int("2"):
62 case str2int("geometric"): {
63 this->printMsg(
64 "Solving with the approximate Dionysus geometric approach.");
65 this->printErr("Not supported");
66 } break;
67 case str2int("3"):
68 case str2int("parallel"): {
69 this->printMsg("Solving with the parallel TTK approach");
70 this->printErr("Not supported");
71 } break;
72 case str2int("bench"): {
73 this->printMsg("Benchmarking");
74 this->printErr("Not supported");
75 } break;
76 default: {
77 this->printErr("You must specify a valid assignment algorithm.");
78 }
79 }
80 }
81
82 this->printMsg("Complete", 1, t.getElapsedTime(), threadNumber_);
83
84 return 0;
85}
86
87void ttk::BottleneckDistance::buildCostMatrices(
88 const ttk::DiagramType &CTDiagram1,
89 const ttk::DiagramType &CTDiagram2,
90 const double zeroThresh,
91 std::vector<std::vector<double>> &minMatrix,
92 std::vector<std::vector<double>> &maxMatrix,
93 std::vector<std::vector<double>> &sadMatrix,
94 const bool reverseMin,
95 const bool reverseMax,
96 const bool reverseSad,
97 const int wasserstein) const {
98
99 int maxI = 0, minI = 0;
100 int maxJ = 0, minJ = 0;
101 int sadI = 0, sadJ = 0;
102
103 for(const auto &p1 : CTDiagram1) {
104 if(std::abs(p1.persistence()) < zeroThresh)
105 continue;
106
107 bool isMin1 = (p1.birth.type == CriticalType::Local_minimum
108 || p1.death.type == CriticalType::Local_minimum);
109 bool isMax1 = (p1.birth.type == CriticalType::Local_maximum
110 || p1.death.type == CriticalType::Local_maximum);
111 bool isSad1 = (p1.birth.type == CriticalType::Saddle1
112 && p1.death.type == CriticalType::Saddle2)
113 || (p1.birth.type == CriticalType::Saddle2
114 && p1.death.type == CriticalType::Saddle1);
115 if(p1.birth.type == CriticalType::Local_minimum
116 && p1.death.type == CriticalType::Local_maximum) {
117 isMin1 = false;
118 isMax1 = true;
119 }
120
121 minJ = 0;
122 maxJ = 0;
123 sadJ = 0;
124
125 for(const auto &p2 : CTDiagram2) {
126 if(std::abs(p2.persistence()) < zeroThresh)
127 continue;
128
129 bool isMin2 = (p2.birth.type == CriticalType::Local_minimum
130 || p2.death.type == CriticalType::Local_minimum);
131 bool isMax2 = (p2.birth.type == CriticalType::Local_maximum
132 || p2.death.type == CriticalType::Local_maximum);
133 bool isSad2 = (p2.birth.type == CriticalType::Saddle1
134 && p2.death.type == CriticalType::Saddle2)
135 || (p2.birth.type == CriticalType::Saddle2
136 && p2.death.type == CriticalType::Saddle1);
137 if(p2.birth.type == CriticalType::Local_minimum
138 && p2.death.type == CriticalType::Local_maximum) {
139 isMin2 = false;
140 isMax2 = true;
141 }
142 if((isMin1 && !isMin2) || (isMax1 && !isMax2) || (isSad1 && !isSad2))
143 continue;
144
145 double distance = this->distanceFunction(p1, p2, wasserstein);
146 double diag1 = this->diagonalDistanceFunction(p1, wasserstein);
147 double diag2 = this->diagonalDistanceFunction(p2, wasserstein);
148
149 if(distance > diag1 + diag2)
150 distance = std::numeric_limits<double>::max();
151
152 if(isMin1 && isMin2) {
153 if(reverseMin)
154 minMatrix[minJ++][minI] = distance;
155 else
156 minMatrix[minI][minJ++] = distance;
157 } else if(isMax1 && isMax2) {
158 if(reverseMax)
159 maxMatrix[maxJ++][maxI] = distance;
160 else
161 maxMatrix[maxI][maxJ++] = distance;
162 } else if(isSad1 && isSad2) {
163 if(reverseSad)
164 sadMatrix[sadJ++][sadI] = distance;
165 else
166 sadMatrix[sadI][sadJ++] = distance;
167 }
168 }
169
170 double distanceToDiagonal = this->diagonalDistanceFunction(p1, wasserstein);
171 if(isMin1) {
172 if(reverseMin)
173 minMatrix[minJ++][minI] = distanceToDiagonal;
174 else
175 minMatrix[minI][minJ++] = distanceToDiagonal;
176 }
177 if(isMax1) {
178 if(reverseMax)
179 maxMatrix[maxJ++][maxI] = distanceToDiagonal;
180 else
181 maxMatrix[maxI][maxJ++] = distanceToDiagonal;
182 }
183 if(isSad1) {
184 if(reverseSad)
185 sadMatrix[sadJ++][sadI] = distanceToDiagonal;
186 else
187 sadMatrix[sadI][sadJ++] = distanceToDiagonal;
188 }
189
190 if(isMin1)
191 ++minI;
192 if(isMax1)
193 ++maxI;
194 if(isSad1)
195 ++sadI;
196 }
197
198 minJ = 0;
199 maxJ = 0;
200 sadJ = 0;
201
202 // Last row: match remaining J components with diagonal.
203 for(const auto &p3 : CTDiagram2) {
204 if(std::abs(p3.persistence()) < zeroThresh)
205 continue;
206
207 bool isMin2 = (p3.birth.type == CriticalType::Local_minimum
208 || p3.death.type == CriticalType::Local_minimum);
209 bool isMax2 = (p3.birth.type == CriticalType::Local_maximum
210 || p3.death.type == CriticalType::Local_maximum);
211 bool isSad2 = (p3.birth.type == CriticalType::Saddle1
212 && p3.death.type == CriticalType::Saddle2)
213 || (p3.birth.type == CriticalType::Saddle2
214 && p3.death.type == CriticalType::Saddle1);
215 if(p3.birth.type == CriticalType::Local_minimum
216 && p3.death.type == CriticalType::Local_maximum) {
217 isMin2 = false;
218 isMax2 = true;
219 }
220
221 double distanceToDiagonal = this->diagonalDistanceFunction(p3, wasserstein);
222 if(isMin2) {
223 if(reverseMin)
224 minMatrix[minJ++][minI] = distanceToDiagonal;
225 else
226 minMatrix[minI][minJ++] = distanceToDiagonal;
227 }
228 if(isMax2) {
229 if(reverseMax)
230 maxMatrix[maxJ++][maxI] = distanceToDiagonal;
231 else
232 maxMatrix[maxI][maxJ++] = distanceToDiagonal;
233 }
234 if(isSad2) {
235 if(reverseSad)
236 sadMatrix[sadJ++][sadI] = distanceToDiagonal;
237 else
238 sadMatrix[sadI][sadJ++] = distanceToDiagonal;
239 }
240 }
241
242 // Last cell
243 if(reverseMin)
244 minMatrix[minJ][minI] = std::numeric_limits<double>::max();
245 else
246 minMatrix[minI][minJ] = std::numeric_limits<double>::max();
247
248 if(reverseMax)
249 maxMatrix[maxJ][maxI] = std::numeric_limits<double>::max();
250 else
251 maxMatrix[maxI][maxJ] = std::numeric_limits<double>::max();
252
253 if(reverseSad)
254 sadMatrix[sadJ][sadI] = std::numeric_limits<double>::max();
255 else
256 sadMatrix[sadI][sadJ] = std::numeric_limits<double>::max();
257}
258
259double ttk::BottleneckDistance::computeGeometricalRange(
260 const ttk::DiagramType &CTDiagram1,
261 const ttk::DiagramType &CTDiagram2) const {
262
263 float minX, minY, minZ, maxX, maxY, maxZ;
264
265 std::array<float, 3> min1{
266 std::numeric_limits<float>::max(),
267 std::numeric_limits<float>::max(),
268 std::numeric_limits<float>::max(),
269 },
270 min2 = min1;
271 std::array<float, 3> max1{
272 std::numeric_limits<float>::lowest(),
273 std::numeric_limits<float>::lowest(),
274 std::numeric_limits<float>::lowest(),
275 },
276 max2 = max1;
277
278 for(const auto &p : CTDiagram1) {
279 min1[0] = std::min(std::min(min1[0], p.birth.coords[0]), p.death.coords[0]);
280 min1[1] = std::min(std::min(min1[1], p.birth.coords[1]), p.death.coords[1]);
281 min1[2] = std::min(std::min(min1[2], p.birth.coords[2]), p.death.coords[2]);
282 max1[0] = std::max(std::max(max1[0], p.birth.coords[0]), p.birth.coords[0]);
283 max1[1] = std::max(std::max(max1[1], p.birth.coords[1]), p.birth.coords[1]);
284 max1[2] = std::max(std::max(max1[2], p.birth.coords[2]), p.birth.coords[2]);
285 }
286
287 for(const auto &p : CTDiagram2) {
288 min2[0] = std::min(std::min(min2[0], p.birth.coords[0]), p.death.coords[0]);
289 min2[1] = std::min(std::min(min2[1], p.birth.coords[1]), p.death.coords[1]);
290 min2[2] = std::min(std::min(min2[2], p.birth.coords[2]), p.death.coords[2]);
291 max2[0] = std::max(std::max(max2[0], p.birth.coords[0]), p.birth.coords[0]);
292 max2[1] = std::max(std::max(max2[1], p.birth.coords[1]), p.birth.coords[1]);
293 max2[2] = std::max(std::max(max2[2], p.birth.coords[2]), p.birth.coords[2]);
294 }
295
296 minX = std::min(min1[0], min2[0]);
297 minY = std::min(min1[1], min2[1]);
298 minZ = std::min(min1[2], min2[2]);
299 maxX = std::max(max1[0], max2[0]);
300 maxY = std::max(max1[1], max2[1]);
301 maxZ = std::max(max1[2], max2[2]);
302
303 const auto square = [](const double a) { return a * a; };
304
305 return std::sqrt(square(maxX - minX) + square(maxY - minY)
306 + square(maxZ - minZ));
307}
308
309double ttk::BottleneckDistance::computeMinimumRelevantPersistence(
310 const ttk::DiagramType &CTDiagram1,
311 const ttk::DiagramType &CTDiagram2) const {
312
313 const auto sp = this->Tolerance;
314 double s = sp > 0.0 && sp < 100.0 ? sp / 100.0 : 0;
315
316 std::vector<double> toSort(CTDiagram1.size() + CTDiagram2.size());
317 for(size_t i = 0; i < CTDiagram1.size(); ++i) {
318 const auto &t = CTDiagram1[i];
319 toSort[i] = std::abs(t.persistence());
320 }
321 for(size_t i = 0; i < CTDiagram2.size(); ++i) {
322 const auto &t = CTDiagram2[i];
323 toSort[CTDiagram1.size() + i] = std::abs(t.persistence());
324 }
325
326 const auto minVal = *std::min_element(toSort.begin(), toSort.end());
327 const auto maxVal = *std::max_element(toSort.begin(), toSort.end());
328 return s * (maxVal - minVal);
329}
330
331void ttk::BottleneckDistance::computeMinMaxSaddleNumberAndMapping(
332 const ttk::DiagramType &CTDiagram,
333 int &nbMin,
334 int &nbMax,
335 int &nbSaddle,
336 std::vector<int> &minMap,
337 std::vector<int> &maxMap,
338 std::vector<int> &sadMap,
339 const double zeroThresh) const {
340
341 for(size_t i = 0; i < CTDiagram.size(); ++i) {
342 const auto &t = CTDiagram[i];
343 const auto nt1 = t.birth.type;
344 const auto nt2 = t.death.type;
345 const auto dt = t.persistence();
346 if(std::abs(dt) < zeroThresh)
347 continue;
348
350 && nt2 == CriticalType::Local_maximum) {
351 nbMax++;
352 maxMap.push_back(i);
353 } else {
355 || nt2 == CriticalType::Local_maximum) {
356 nbMax++;
357 maxMap.push_back(i);
358 }
360 || nt2 == CriticalType::Local_minimum) {
361 nbMin++;
362 minMap.push_back(i);
363 }
364 if((nt1 == CriticalType::Saddle1 && nt2 == CriticalType::Saddle2)
365 || (nt1 == CriticalType::Saddle2 && nt2 == CriticalType::Saddle1)) {
366 nbSaddle++;
367 sadMap.push_back(i);
368 }
369 }
370 }
371}
372
373static void solvePWasserstein(std::vector<std::vector<double>> &matrix,
374 std::vector<ttk::MatchingType> &matchings,
376
377 solver.setInput(matrix);
378 solver.run(matchings);
379 solver.clearMatrix();
380}
381
382static void solveInfinityWasserstein(const int nbRow,
383 const int nbCol,
384 std::vector<std::vector<double>> &matrix,
385 std::vector<ttk::MatchingType> &matchings,
386 ttk::GabowTarjan &solver) {
387
388 // Copy input matrix.
389 auto bottleneckMatrix = matrix;
390
391 // Solve.
392 solver.setInput(nbRow, nbCol, &bottleneckMatrix);
393 solver.run(matchings);
394 solver.clear();
395}
396
397double ttk::BottleneckDistance::buildMappings(
398 const std::vector<MatchingType> &inputMatchings,
399 const bool transposeGlobal,
400 const bool transposeLocal,
401 std::vector<MatchingType> &outputMatchings,
402 const std::vector<int> &m1,
403 const std::vector<int> &m2,
404 int wasserstein) const {
405
406 // Input map permutation (so as to ignore transposition later on)
407 const auto &map1 = transposeLocal ? m2 : m1;
408 const auto &map2 = transposeLocal ? m1 : m2;
409
410 double addedPersistence = 0;
411 const auto doTranspose = transposeGlobal ^ transposeLocal;
412
413 for(const auto &t : inputMatchings) {
414
415 const auto val = std::abs(std::get<2>(t));
416 const size_t p1 = std::get<0>(t);
417 const size_t p2 = std::get<1>(t);
418
419 if(p1 >= map1.size() || p2 >= map2.size()) {
420 addedPersistence = (wasserstein > 0 ? addedPersistence + val
421 : std::max(val, addedPersistence));
422 } else {
423 if(doTranspose) {
424 outputMatchings.emplace_back(map2[p2], map1[p1], val);
425 } else {
426 outputMatchings.emplace_back(map1[p1], map2[p2], val);
427 }
428 }
429 }
430
431 return addedPersistence;
432}
433
434double ttk::BottleneckDistance::distanceFunction(const ttk::PersistencePair &a,
435 const ttk::PersistencePair &b,
436 const int wasserstein) const {
437
438 const int w = std::max(wasserstein, 1); // L_inf not managed.
439
440 // We don't match critical points of different index.
441 // This must be ensured before calling the distance function.
442 const bool isMin1 = a.birth.type == CriticalType::Local_minimum;
443 const bool isMax1 = a.death.type == CriticalType::Local_maximum;
444
445 const std::array<float, 3> coordsAbsDiff{
446 std::abs(a.death.coords[0] - b.death.coords[0]),
447 std::abs(a.death.coords[1] - b.death.coords[1]),
448 std::abs(a.death.coords[2] - b.death.coords[2]),
449 };
450
451 const auto x
452 = ((isMin1 && !isMax1) ? this->PE : this->PS)
453 * Geometry::pow(std::abs(a.birth.sfValue - b.birth.sfValue), w);
454 const auto y
455 = (isMax1 ? this->PE : this->PS)
456 * Geometry::pow(std::abs(a.death.sfValue - b.death.sfValue), w);
457 const double geoDistance
458 = isMax1 || isMin1
459 ? (this->PX * Geometry::pow(coordsAbsDiff[0], w)
460 + this->PY * Geometry::pow(coordsAbsDiff[1], w)
461 + this->PZ * Geometry::pow(coordsAbsDiff[2], w))
462 : (this->PX
464 std::abs(a.birth.coords[0] + a.death.coords[0]) / 2
465 - std::abs(b.birth.coords[0] + b.death.coords[0]) / 2,
466 w)
467 + this->PY
469 std::abs(a.birth.coords[1] + a.death.coords[1]) / 2
470 - std::abs(b.birth.coords[1] + b.death.coords[1]) / 2,
471 w)
472 + this->PZ
474 std::abs(a.birth.coords[2] + a.death.coords[2]) / 2
475 - std::abs(b.birth.coords[2] + b.death.coords[2]) / 2,
476 w));
477
478 double persDistance = x + y;
479 return Geometry::pow(persDistance + geoDistance, 1.0 / w);
480}
481
482double ttk::BottleneckDistance::diagonalDistanceFunction(
483 const ttk::PersistencePair &a, const int wasserstein) const {
484
485 const int w = std::max(wasserstein, 1);
486 const bool isMin1 = a.birth.type == CriticalType::Local_minimum;
487 const bool isMax1 = a.death.type == CriticalType::Local_maximum;
488 const double infDistance = (isMin1 || isMax1 ? this->PE : this->PS)
489 * Geometry::pow(std::abs(a.persistence()), w);
490 const double geoDistance
491 = (this->PX
492 * Geometry::pow(std::abs(a.death.coords[0] - a.birth.coords[0]), w)
493 + this->PY
494 * Geometry::pow(std::abs(a.death.coords[1] - a.birth.coords[1]), w)
495 + this->PZ
496 * Geometry::pow(std::abs(a.death.coords[2] - a.birth.coords[2]), w));
497
498 return Geometry::pow(infDistance + geoDistance, 1.0 / w);
499}
500
501int ttk::BottleneckDistance::computeBottleneck(
502 const ttk::DiagramType &d1,
503 const ttk::DiagramType &d2,
504 std::vector<MatchingType> &matchings) {
505
506 const auto transposeOriginal = d1.size() > d2.size();
507 if(transposeOriginal) {
508 this->printMsg("The first persistence diagram is larger than the second.");
509 this->printMsg("Solving the transposed problem.");
510 }
511 const auto &CTDiagram1 = transposeOriginal ? d2 : d1;
512 const auto &CTDiagram2 = transposeOriginal ? d1 : d2;
513
514 // Check user parameters.
515 const auto isBottleneck = this->WassersteinMetric == "inf";
516 const int wasserstein = isBottleneck ? -1 : stoi(this->WassersteinMetric);
517 if(wasserstein < 0 && !isBottleneck) {
518 this->printErr("Wrong value for the Wassertein power parameter");
519 return -4;
520 }
521
522 // Needed to limit computation time.
523 const auto zeroThresh
524 = this->computeMinimumRelevantPersistence(CTDiagram1, CTDiagram2);
525
526 // Initialize solvers.
527 std::vector<MatchingType> minMatchings;
528 std::vector<MatchingType> maxMatchings;
529 std::vector<MatchingType> sadMatchings;
530
531 // Initialize cost matrices.
532 int nbRowMin = 0, nbColMin = 0;
533 int nbRowMax = 0, nbColMax = 0;
534 int nbRowSad = 0, nbColSad = 0;
535
536 // Remap for matchings.
537 std::array<std::vector<int>, 3> map1{};
538 std::array<std::vector<int>, 3> map2{};
539
540 this->computeMinMaxSaddleNumberAndMapping(CTDiagram1, nbRowMin, nbRowMax,
541 nbRowSad, map1[0], map1[2], map1[1],
542 zeroThresh);
543 this->computeMinMaxSaddleNumberAndMapping(CTDiagram2, nbColMin, nbColMax,
544 nbColSad, map2[0], map2[2], map2[1],
545 zeroThresh);
546
547 // Automatically transpose if nb rows > nb cols
548 const auto maxRowColMin = std::max(nbRowMin + 1, nbColMin + 1);
549 const auto maxRowColMax = std::max(nbRowMax + 1, nbColMax + 1);
550 const auto maxRowColSad = std::max(nbRowSad + 1, nbColSad + 1);
551
552 const auto minRowColMin = std::min(nbRowMin + 1, nbColMin + 1);
553 const auto minRowColMax = std::min(nbRowMax + 1, nbColMax + 1);
554 const auto minRowColSad = std::min(nbRowSad + 1, nbColSad + 1);
555
556 std::vector<std::vector<double>> minMatrix(
557 minRowColMin, std::vector<double>(maxRowColMin));
558 std::vector<std::vector<double>> maxMatrix(
559 minRowColMax, std::vector<double>(maxRowColMax));
560 std::vector<std::vector<double>> sadMatrix(
561 minRowColSad, std::vector<double>(maxRowColSad));
562
563 const bool transposeMin = nbRowMin > nbColMin;
564 const bool transposeMax = nbRowMax > nbColMax;
565 const bool transposeSad = nbRowSad > nbColSad;
566
567 Timer t;
568
569 this->buildCostMatrices(CTDiagram1, CTDiagram2, zeroThresh, minMatrix,
570 maxMatrix, sadMatrix, transposeMin, transposeMax,
571 transposeSad, wasserstein);
572
573 if(!isBottleneck) {
574
575 if(nbRowMin > 0 && nbColMin > 0) {
576 AssignmentMunkres<double> solverMin;
577 this->printMsg("Affecting minima...");
578 solvePWasserstein(minMatrix, minMatchings, solverMin);
579 }
580
581 if(nbRowMax > 0 && nbColMax > 0) {
582 AssignmentMunkres<double> solverMax;
583 this->printMsg("Affecting maxima...");
584 solvePWasserstein(maxMatrix, maxMatchings, solverMax);
585 }
586
587 if(nbRowSad > 0 && nbColSad > 0) {
588 AssignmentMunkres<double> solverSad;
589 this->printMsg("Affecting saddles...");
590 solvePWasserstein(sadMatrix, sadMatchings, solverSad);
591 }
592
593 } else {
594
595 // Launch solving for minima.
596 if(nbRowMin > 0 && nbColMin > 0) {
597 GabowTarjan solverMin;
598 this->printMsg("Affecting minima...");
599 solveInfinityWasserstein(
600 minRowColMin, maxRowColMin, minMatrix, minMatchings, solverMin);
601 }
602
603 // Launch solving for maxima.
604 if(nbRowMax > 0 && nbColMax > 0) {
605 GabowTarjan solverMax;
606 this->printMsg("Affecting maxima...");
607 solveInfinityWasserstein(
608 minRowColMax, maxRowColMax, maxMatrix, maxMatchings, solverMax);
609 }
610
611 // Launch solving for saddles.
612 if(nbRowSad > 0 && nbColSad > 0) {
613 GabowTarjan solverSad;
614 this->printMsg("Affecting saddles...");
615 solveInfinityWasserstein(
616 minRowColSad, maxRowColSad, sadMatrix, sadMatchings, solverSad);
617 }
618 }
619
620 this->printMsg("TTK CORE DONE", 1, t.getElapsedTime());
621
622 // Rebuild mappings.
623 // Begin cost computation for unpaired vertices.
624 const std::array<double, 3> addedPersistence{
625 this->buildMappings(minMatchings, transposeOriginal, transposeMin,
626 matchings, map1[0], map2[0], wasserstein),
627 this->buildMappings(sadMatchings, transposeOriginal, transposeSad,
628 matchings, map1[1], map2[1], wasserstein),
629 this->buildMappings(maxMatchings, transposeOriginal, transposeMax,
630 matchings, map1[2], map2[2], wasserstein),
631 };
632
633 // TODO [HIGH] do that for embeddings
634 // Recompute matching weights for user-friendly distance.
635 std::array<double, 3> costs{};
636 std::vector<bool> paired1(CTDiagram1.size(), false);
637 std::vector<bool> paired2(CTDiagram2.size(), false);
638
639 // store last matching distance for bottleneck
640 double partialDistance{};
641
642 for(const auto &mt : matchings) {
643 int i = transposeOriginal ? std::get<1>(mt) : std::get<0>(mt);
644 int j = transposeOriginal ? std::get<0>(mt) : std::get<1>(mt);
645
646 const auto &t1 = CTDiagram1[i];
647 const auto &t2 = CTDiagram2[j];
648 paired1[i] = true;
649 paired2[j] = true;
650
651 partialDistance = this->distanceFunction(t1, t2, wasserstein);
652
653 if(t1.death.type == CriticalType::Local_maximum) {
654 if(!isBottleneck) {
655 costs[2] += partialDistance;
656 } else {
657 costs[2] = std::max(costs[2], partialDistance);
658 }
659 } else if(t1.birth.type == CriticalType::Local_minimum) {
660 if(!isBottleneck) {
661 costs[0] += partialDistance;
662 } else {
663 costs[0] = std::max(costs[0], partialDistance);
664 }
665 } else if(t1.birth.type == CriticalType::Saddle1
666 && t1.death.type == CriticalType::Saddle2) {
667 if(!isBottleneck) {
668 costs[1] += partialDistance;
669 } else {
670 costs[1] = std::max(costs[1], partialDistance);
671 }
672 }
673 }
674
675 const auto affectationD
676 = !isBottleneck ? costs[0] + costs[1] + costs[2] : partialDistance;
677 const auto addedPers
678 = addedPersistence[0] + addedPersistence[1] + addedPersistence[2];
679 this->distance_
680 = !isBottleneck
681 ? Geometry::pow(affectationD + addedPers, 1.0 / wasserstein)
682 : std::max(affectationD, *std::max_element(addedPersistence.begin(),
683 addedPersistence.end()));
684
685 std::stringstream msg;
686 this->printMsg("Computed distance:");
687 this->printMsg("diagMax(" + std::to_string(addedPersistence[2])
688 + "), diagMin(" + std::to_string(addedPersistence[0])
689 + "), diagSad(" + std::to_string(addedPersistence[1]) + ")");
690 this->printMsg("affAll(" + std::to_string(affectationD) + "), res("
691 + std::to_string(this->distance_) + ")");
692
693 // aggregate costs per pair type
694 if(!isBottleneck) {
695 costs[0] = Geometry::pow(costs[0] + addedPersistence[0], 1.0 / wasserstein);
696 costs[1] = Geometry::pow(costs[1] + addedPersistence[1], 1.0 / wasserstein);
697 costs[2] = Geometry::pow(costs[2] + addedPersistence[2], 1.0 / wasserstein);
698 } else {
699 costs[0] = std::max(costs[0], addedPersistence[0]);
700 costs[1] = std::max(costs[1], addedPersistence[1]);
701 costs[2] = std::max(costs[2], addedPersistence[2]);
702 }
703
704 this->costs_ = costs;
705
706 // display results
707 std::vector<std::vector<std::string>> rows{
708 {" Min-saddle cost", std::to_string(this->costs_[0])},
709 {" Saddle-saddle cost", std::to_string(this->costs_[1])},
710 {" Saddle-max cost", std::to_string(this->costs_[2])},
711 {isBottleneck ? "Bottleneck Distance" : "Wasserstein Distance",
712 std::to_string(this->distance_)},
713 };
714 this->printMsg(rows);
715
716 return 0;
717}
int run(std::vector< MatchingType > &matchings) override
int setInput(std::vector< std::vector< dataType > > &C_) override
virtual void clearMatrix()
int execute(const ttk::DiagramType &diag0, const ttk::DiagramType &diag1, std::vector< MatchingType > &matchings)
void setDebugMsgPrefix(const std::string &prefix)
Definition: Debug.h:364
int run(std::vector< MatchingType > &matchings)
void setInput(int rowSize_, int colSize_, std::vector< std::vector< double > > *C_)
Definition: GabowTarjan.h:40
double getElapsedTime()
Definition: Timer.h:15
T1 pow(const T1 val, const T2 n)
Definition: Geometry.h:411
T distance(const T *p0, const T *p1, const int &dimension=3)
Definition: Geometry.cpp:344
constexpr unsigned long long str2int(const char *str, int h=0)
std::vector< PersistencePair > DiagramType
Persistence Diagram type as a vector of Persistence pairs.
std::array< float, 3 > coords
ttk::CriticalVertex birth
double persistence() const
Return the topological persistence of the pair.
ttk::CriticalVertex death
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)