TTK
Loading...
Searching...
No Matches
PersistenceDiagramDistanceMatrix.cpp
Go to the documentation of this file.
1#include <algorithm>
2#include <limits>
3
5
6using namespace ttk;
7
8std::vector<std::vector<double>> PersistenceDiagramDistanceMatrix::execute(
9 const std::vector<DiagramType> &intermediateDiagrams,
10 const std::array<size_t, 2> &nInputs) const {
11
12 Timer tm{};
13
14 const auto nDiags = intermediateDiagrams.size();
15
16 if(do_min_ && do_sad_ && do_max_) {
17 this->printMsg("Processing all critical pairs types");
18 } else if(do_min_) {
19 this->printMsg("Processing only MIN-SAD pairs");
20 } else if(do_sad_) {
21 this->printMsg("Processing only SAD-SAD pairs");
22 } else if(do_max_) {
23 this->printMsg("Processing only SAD-MAX pairs");
24 }
25
26 std::vector<DiagramType> inputDiagramsMin(nDiags);
27 std::vector<DiagramType> inputDiagramsSad(nDiags);
28 std::vector<DiagramType> inputDiagramsMax(nDiags);
29
30 std::vector<BidderDiagram> bidder_diagrams_min{};
31 std::vector<BidderDiagram> bidder_diagrams_sad{};
32 std::vector<BidderDiagram> bidder_diagrams_max{};
33 std::vector<BidderDiagram> current_bidder_diagrams_min{};
34 std::vector<BidderDiagram> current_bidder_diagrams_sad{};
35 std::vector<BidderDiagram> current_bidder_diagrams_max{};
36
37 // Store the persistence of the global min-max pair
38 std::vector<double> maxDiagPersistence(nDiags);
39
40 // Create diagrams for min, saddle and max persistence pairs
41#ifdef TTK_ENABLE_OPENMP
42#pragma omp parallel for num_threads(threadNumber_)
43#endif // TTK_ENABLE_OPENMP
44 for(size_t i = 0; i < nDiags; i++) {
45 for(const auto &p : intermediateDiagrams[i]) {
46 maxDiagPersistence[i] = std::max(p.persistence(), maxDiagPersistence[i]);
47
48 if(p.persistence() > 0) {
49 if(p.birth.type == CriticalType::Local_minimum
50 && p.death.type == CriticalType::Local_maximum) {
51 inputDiagramsMax[i].emplace_back(p);
52 } else {
53 if(p.birth.type == CriticalType::Local_maximum
54 || p.death.type == CriticalType::Local_maximum) {
55 inputDiagramsMax[i].emplace_back(p);
56 }
57 if(p.birth.type == CriticalType::Local_minimum
58 || p.death.type == CriticalType::Local_minimum) {
59 inputDiagramsMin[i].emplace_back(p);
60 }
61 if((p.birth.type == CriticalType::Saddle1
62 && p.death.type == CriticalType::Saddle2)
63 || (p.birth.type == CriticalType::Saddle2
64 && p.death.type == CriticalType::Saddle1)) {
65 inputDiagramsSad[i].emplace_back(p);
66 }
67 }
68 }
69 }
70 }
71
72 if(this->do_min_) {
73 setBidderDiagrams(nDiags, inputDiagramsMin, bidder_diagrams_min);
74 }
75 if(this->do_sad_) {
76 setBidderDiagrams(nDiags, inputDiagramsSad, bidder_diagrams_sad);
77 }
78 if(this->do_max_) {
79 setBidderDiagrams(nDiags, inputDiagramsMax, bidder_diagrams_max);
80 }
81
82 switch(this->Constraint) {
84 this->printMsg("Using all diagram pairs");
85 break;
87 this->printMsg("Using the " + std::to_string(this->MaxNumberOfPairs)
88 + " most persistent pairs");
89 break;
91 std::stringstream pers{};
92 pers << std::fixed << std::setprecision(2) << this->MinPersistence;
93 this->printMsg("Using diagram pairs above a persistence threshold of "
94 + pers.str());
95 } break;
97 this->printMsg(
98 "Using the "
99 + std::to_string(static_cast<int>(100 * (1 - this->MinPersistence)))
100 + "% most persistent pairs of every diagram");
101 break;
103 this->printMsg(
104 "Using the "
105 + std::to_string(static_cast<int>(100 * (1 - this->MinPersistence)))
106 + "% most persistent pairs of all diagrams");
107 break;
108 }
109
110 std::vector<std::vector<double>> distMat{};
112 getDiagramsDistMat(nInputs, distMat, bidder_diagrams_min,
113 bidder_diagrams_sad, bidder_diagrams_max);
114 } else {
115 if(this->do_min_) {
117 bidder_diagrams_min, current_bidder_diagrams_min, maxDiagPersistence);
118 }
119 if(this->do_sad_) {
121 bidder_diagrams_sad, current_bidder_diagrams_sad, maxDiagPersistence);
122 }
123 if(this->do_max_) {
125 bidder_diagrams_max, current_bidder_diagrams_max, maxDiagPersistence);
126 }
127 getDiagramsDistMat(nInputs, distMat, current_bidder_diagrams_min,
128 current_bidder_diagrams_sad,
129 current_bidder_diagrams_max);
130 }
131
132 this->printMsg("Complete", 1.0, tm.getElapsedTime(), this->threadNumber_);
133
134 return distMat;
135}
136
138 const std::vector<BidderDiagram> &bidder_diags) const {
139
140 double max_persistence = 0;
141
142 for(unsigned int i = 0; i < bidder_diags.size(); ++i) {
143 for(size_t j = 0; j < bidder_diags[i].size(); ++j) {
144 const double persistence = bidder_diags[i][j].getPersistence();
145 if(persistence > max_persistence) {
146 max_persistence = persistence;
147 }
148 }
149 }
150
151 return max_persistence;
152}
153
155 const BidderDiagram &D1, const BidderDiagram &D2) const {
156
157 GoodDiagram D2_bis{};
158 for(size_t i = 0; i < D2.size(); i++) {
159 const Bidder &b = D2[i];
160 Good g(b.x_, b.y_, b.isDiagonal(), D2_bis.size());
162 g.setPrice(0);
163 D2_bis.emplace_back(g);
164 }
165
167 this->Wasserstein, this->Alpha, this->Lambda, this->DeltaLim, true);
168 auction.BuildAuctionDiagrams(D1, D2_bis);
169 return auction.run();
170}
171
173 const std::array<size_t, 2> &nInputs,
174 std::vector<std::vector<double>> &distanceMatrix,
175 const std::vector<BidderDiagram> &diags_min,
176 const std::vector<BidderDiagram> &diags_sad,
177 const std::vector<BidderDiagram> &diags_max) const {
178
179 distanceMatrix.resize(nInputs[0]);
180
181#ifdef TTK_ENABLE_OPENMP
182#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
183#endif // TTK_ENABLE_OPENMP
184 for(size_t i = 0; i < nInputs[0]; ++i) {
185
186 if(nInputs[1] == 0) {
187 distanceMatrix[i].resize(nInputs[0]);
188 // set the matrix diagonal
189 distanceMatrix[i][i] = 0.0;
190 } else {
191 distanceMatrix[i].resize(nInputs[1]);
192 }
193
194 const auto getDist = [&](const size_t a, const size_t b) -> double {
195 double distance{};
196 if(this->do_min_) {
197 auto &dimin = diags_min[a];
198 auto &djmin = diags_min[b];
199 distance += computePowerDistance(dimin, djmin);
200 }
201 if(this->do_sad_) {
202 auto &disad = diags_sad[a];
203 auto &djsad = diags_sad[b];
204 distance += computePowerDistance(disad, djsad);
205 }
206 if(this->do_max_) {
207 auto &dimax = diags_max[a];
208 auto &djmax = diags_max[b];
209 distance += computePowerDistance(dimax, djmax);
210 }
211 return Geometry::pow(distance, 1.0 / this->Wasserstein);
212 };
213
214 if(nInputs[1] == 0) {
215 // square matrix: only compute the upper triangle (i < j < nInputs[0])
216 for(size_t j = i + 1; j < nInputs[0]; ++j) {
217 distanceMatrix[i][j] = getDist(i, j);
218 }
219 } else {
220 // rectangular matrix: compute the whole line/column (0 <= j < nInputs[1])
221 for(size_t j = 0; j < nInputs[1]; ++j) {
222 distanceMatrix[i][j] = getDist(i, j + nInputs[0]);
223 }
224 }
225 }
226
227 if(nInputs[1] == 0) {
228 // square distance matrix is symmetric: complete the lower triangle
229 for(size_t i = 0; i < nInputs[0]; ++i) {
230 for(size_t j = i + 1; j < nInputs[0]; ++j) {
231 distanceMatrix[j][i] = distanceMatrix[i][j];
232 }
233 }
234 }
235}
236
238 const size_t nInputs,
239 std::vector<DiagramType> &inputDiagrams,
240 std::vector<BidderDiagram> &bidder_diags) const {
241
242 bidder_diags.resize(nInputs);
243
244 for(size_t i = 0; i < nInputs; i++) {
245 auto &diag = inputDiagrams[i];
246 auto &bidders = bidder_diags[i];
247
248 for(size_t j = 0; j < diag.size(); j++) {
249 // Add bidder to bidders
250 Bidder b(diag[j], j, this->Lambda);
251 b.setPositionInAuction(bidders.size());
252 bidders.emplace_back(b);
253 if(b.isDiagonal() || b.x_ == b.y_) {
254 this->printMsg("Diagonal point in diagram " + std::to_string(i) + "!",
256 }
257 }
258 }
259}
260
262 const std::vector<BidderDiagram> &bidder_diags,
263 std::vector<BidderDiagram> &current_bidder_diags,
264 const std::vector<double> &maxDiagPersistence) const {
265
266 current_bidder_diags.resize(bidder_diags.size());
267 const auto nInputs = current_bidder_diags.size();
268 const auto maxPersistence
269 = *std::max_element(maxDiagPersistence.begin(), maxDiagPersistence.end());
270
274 for(size_t i = 0; i < nInputs; ++i) {
275 for(auto b : bidder_diags[i]) {
276
277 if( // filter out pairs below absolute persistence threshold
279 && b.getPersistence() > this->MinPersistence)
280 || // filter out pairs below persistence threshold relative to
281 // the most persistent pair *of each diagrams*
283 && b.getPersistence() > this->MinPersistence * maxDiagPersistence[i])
284 || // filter out pairs below persistence threshold relative to the
285 // most persistence pair *in all diagrams*
287 && b.getPersistence() > this->MinPersistence * maxPersistence)) {
288 b.id_ = current_bidder_diags[i].size();
289 b.setPositionInAuction(current_bidder_diags[i].size());
290 current_bidder_diags[i].emplace_back(b);
291 }
292 }
293 }
294 return;
295 }
296
297 const double prev_min_persistence = 2.0 * getMostPersistent(bidder_diags);
298 double new_min_persistence = 0.0;
299
300 // 1. Get size of the largest current diagram, deduce the maximal number
301 // of points to append
302 size_t max_diagram_size = 0;
303 for(const auto &diag : current_bidder_diags) {
304 max_diagram_size = std::max(diag.size(), max_diagram_size);
305 }
306 size_t const max_points_to_add = std::max(
307 this->MaxNumberOfPairs, this->MaxNumberOfPairs + max_diagram_size / 10);
308 // 2. Get which points can be added, deduce the new minimal persistence
309 std::vector<std::vector<int>> candidates_to_be_added(nInputs);
310 std::vector<std::vector<size_t>> idx(nInputs);
311
312 for(size_t i = 0; i < nInputs; i++) {
313 double local_min_persistence = std::numeric_limits<double>::min();
314 std::vector<double> persistences;
315 for(size_t j = 0; j < bidder_diags[i].size(); j++) {
316 const auto &b = bidder_diags[i][j];
317 double const persistence = b.getPersistence();
318 if(persistence >= 0.0 && persistence <= prev_min_persistence) {
319 candidates_to_be_added[i].emplace_back(j);
320 idx[i].emplace_back(idx[i].size());
321 persistences.emplace_back(persistence);
322 }
323 }
324 const auto cmp = [&persistences](const size_t a, const size_t b) {
325 return ((persistences[a] > persistences[b])
326 || ((persistences[a] == persistences[b]) && (a > b)));
327 };
328 std::sort(idx[i].begin(), idx[i].end(), cmp);
329 const auto size = candidates_to_be_added[i].size();
330 if(size >= max_points_to_add) {
331 double const last_persistence_added
332 = persistences[idx[i][max_points_to_add - 1]];
333 if(last_persistence_added > local_min_persistence) {
334 local_min_persistence = last_persistence_added;
335 }
336 }
337 if(i == 0) {
338 new_min_persistence = local_min_persistence;
339 } else {
340 if(local_min_persistence < new_min_persistence) {
341 new_min_persistence = local_min_persistence;
342 }
343 }
344 // 3. Add the points to the current diagrams
345 const auto s = candidates_to_be_added[i].size();
346 for(size_t j = 0; j < std::min(max_points_to_add, s); j++) {
347 auto b = bidder_diags[i].at(candidates_to_be_added[i][idx[i][j]]);
348 const double persistence = b.getPersistence();
349 if(persistence >= new_min_persistence) {
350 b.id_ = current_bidder_diags[i].size();
351 b.setPositionInAuction(current_bidder_diags[i].size());
352 current_bidder_diags[i].emplace_back(b);
353 }
354 }
355 }
356}
void setPositionInAuction(const int pos)
void setPrice(const double price)
void SetCriticalCoordinates(const float coords_x, const float coords_y, const float coords_z)
double run(std::vector< MatchingType > &matchings, const int kdt_index=0)
void BuildAuctionDiagrams(const BidderDiagram &BD, const GoodDiagram &GD)
void enrichCurrentBidderDiagrams(const std::vector< BidderDiagram > &bidder_diags, std::vector< BidderDiagram > &current_bidder_diags, const std::vector< double > &maxDiagPersistence) const
std::vector< std::vector< double > > execute(const std::vector< DiagramType > &intermediateDiagrams, const std::array< size_t, 2 > &nInputs) const
double getMostPersistent(const std::vector< BidderDiagram > &bidder_diags) const
void setBidderDiagrams(const size_t nInputs, std::vector< DiagramType > &inputDiagrams, std::vector< BidderDiagram > &bidder_diags) const
double computePowerDistance(const BidderDiagram &D1, const BidderDiagram &D2) const
void getDiagramsDistMat(const std::array< size_t, 2 > &nInputs, std::vector< std::vector< double > > &distanceMatrix, const std::vector< BidderDiagram > &diags_min, const std::vector< BidderDiagram > &diags_sad, const std::vector< BidderDiagram > &diags_max) const
T1 pow(const T1 val, const T2 n)
Definition Geometry.h:456
The Topology ToolKit.
std::vector< Bidder > BidderDiagram
std::vector< Good > GoodDiagram
T end(std::pair< T, T > &p)
Definition ripserpy.cpp:483
T begin(std::pair< T, T > &p)
Definition ripserpy.cpp:479
printMsg(debug::output::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/|__ _|"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)