TTK
Loading...
Searching...
No Matches
DiscreteMorseSandwich.cpp
Go to the documentation of this file.
2
4 this->setDebugMsgPrefix("DiscreteMorseSandwich");
5}
6
8 std::vector<PersistencePair> &pairs,
9 std::vector<bool> &pairedExtrema,
10 std::vector<bool> &pairedSaddles,
11 std::vector<SimplexId> &reps,
12 std::vector<tripletType> &triplets,
13 const SimplexId *const saddlesOrder,
14 const SimplexId *const extremaOrder,
15 const SimplexId pairDim) const {
16
17 // comparison functions
18 const auto cmpSadMax
19 = [=](const tripletType &t0, const tripletType &t1) -> bool {
20 const auto s0 = t0[0];
21 const auto s1 = t1[0];
22 const auto m0 = t0[2];
23 const auto m1 = t1[2];
24
25#ifdef _LIBCPP_VERSION
26 // libc++'s std::sort compares an entry to itself
27 if(&t0 == &t1) {
28 return true;
29 }
30#endif // _LIBCPP_VERSION
31
32 if(s0 != s1)
33 return saddlesOrder[s0] > saddlesOrder[s1];
34 else
35 return extremaOrder[m0] < extremaOrder[m1];
36 };
37
38 const auto cmpSadMin
39 = [=](const tripletType &t0, const tripletType &t1) -> bool {
40 const auto s0 = t0[0];
41 const auto s1 = t1[0];
42 const auto m0 = t0[2];
43 const auto m1 = t1[2];
44 if(s0 != s1)
45 return saddlesOrder[s0] < saddlesOrder[s1];
46 else
47 return extremaOrder[m0] > extremaOrder[m1];
48 };
49
50 // sort triplets
51 if(pairDim == 0) {
52 TTK_PSORT(this->threadNumber_, triplets.begin(), triplets.end(), cmpSadMin);
53 } else {
54 // saddle-saddle pairs from 1-saddles to 2-saddles
55 TTK_PSORT(this->threadNumber_, triplets.begin(), triplets.end(), cmpSadMax);
56 }
57
58 // get representative of current extremum
59 const auto getRep = [&reps](SimplexId v) -> SimplexId {
60 auto r = reps[v];
61 while(r != v) {
62 v = r;
63 r = reps[v];
64 }
65 return r;
66 };
67
68 const bool increasing = (pairDim > 0);
69
70 const auto addPair = [&pairs, &pairedExtrema, &pairedSaddles, increasing,
71 pairDim](const SimplexId sad, const SimplexId extr) {
72 if(increasing) {
73 pairs.emplace_back(sad, extr, pairDim);
74 } else {
75 pairs.emplace_back(extr, sad, pairDim);
76 }
77 pairedSaddles[sad] = true;
78 pairedExtrema[extr] = true;
79 };
80
81 for(const auto &t : triplets) {
82 const auto sv = t[0];
83
84 auto r1 = getRep(t[1]);
85
86 if(t[2] < 0) {
87 // deal with "shadow" triplets (a 2-saddle with only one
88 // ascending 1-separatrix leading to an unique maximum)
89 if(!pairedExtrema[r1] && !pairedSaddles[sv]) {
90 // when considering the boundary, the "-1" of the triplets
91 // indicate a virtual maximum of infinite persistence on the
92 // boundary component. a pair is created with the other
93 // maximum
94 addPair(sv, r1);
95 }
96
97 continue;
98 }
99
100 auto r2 = getRep(t[2]);
101
102 if(r1 != r2) {
103 if(((extremaOrder[r1] > extremaOrder[r2]) == increasing
104 || pairedExtrema[r1])
105 && !pairedExtrema[r2]) {
106 std::swap(r1, r2);
107 }
108 if(!pairedExtrema[r1]) {
109 addPair(sv, r1);
110 reps[t[1]] = r2;
111 reps[r1] = r2;
112 }
113 }
114 }
115}
116
118 const std::vector<PersistencePair> &pairs,
119 const std::array<std::vector<SimplexId>, 4> &criticalCellsByDim,
120 const std::vector<bool> &pairedMinima,
121 const std::vector<bool> &paired1Saddles,
122 const std::vector<bool> &paired2Saddles,
123 const std::vector<bool> &pairedMaxima) const {
124
125 const auto dim = this->dg_.getDimensionality();
126
127 // display number of pairs per pair type
128 std::vector<std::vector<std::string>> rows{
129 {" #Min-saddle pairs",
130 std::to_string(
131 std::count_if(pairs.begin(), pairs.end(),
132 [](const PersistencePair &a) { return a.type == 0; }))},
133 {" #Saddle-saddle pairs",
134 std::to_string(dim == 3 ? std::count_if(
135 pairs.begin(), pairs.end(),
136 [](const PersistencePair &a) { return a.type == 1; })
137 : 0)},
138 {" #Saddle-max pairs",
139 std::to_string(std::count_if(
140 pairs.begin(), pairs.end(),
141 [dim](const PersistencePair &a) { return a.type == dim - 1; }))},
142 };
143
144 // display number of critical cells (paired and unpaired)
145 std::vector<size_t> nCritCells(dim + 1);
146 std::vector<size_t> nNonPairedCritCells(dim + 1);
147
148 for(int i = 0; i < dim + 1; ++i) {
149 nCritCells[i] = criticalCellsByDim[i].size();
150 size_t nNonPaired{};
151 for(size_t j = 0; j < criticalCellsByDim[i].size(); ++j) {
152 const auto cell = criticalCellsByDim[i][j];
153 if((i == 0 && !pairedMinima[cell]) || (i == 1 && !paired1Saddles[cell])
154 || (i == 2 && dim == 3 && !paired2Saddles[cell])
155 || (i == dim && !pairedMaxima[cell])) {
156 nNonPaired++;
157 }
158 }
159 nNonPairedCritCells[i] = nNonPaired;
160 }
161
162 std::vector<std::string> critCellsLabels{"Minima"};
163 if(dim >= 2) {
164 critCellsLabels.emplace_back("1-saddles");
165 }
166 if(dim >= 3) {
167 critCellsLabels.emplace_back("2-saddles");
168 }
169 critCellsLabels.emplace_back("Maxima");
170
171 for(int i = 0; i < dim + 1; ++i) {
172 const std::string unpaired{nNonPairedCritCells[i] == 0
173 ? " (all paired)"
174 : " (" + std::to_string(nNonPairedCritCells[i])
175 + " unpaired)"};
176
177 rows.emplace_back(std::vector<std::string>{
178 " #" + critCellsLabels[i], std::to_string(nCritCells[i]) + unpaired});
179 }
181}
#define TTK_PSORT(NTHREADS,...)
Parallel sort macro.
Definition: OpenMP.h:46
void setDebugMsgPrefix(const std::string &prefix)
Definition: Debug.h:364
std::array< SimplexId, 3 > tripletType
Triplet type for persistence pairs.
void displayStats(const std::vector< PersistencePair > &pairs, const std::array< std::vector< SimplexId >, 4 > &criticalCellsByDim, const std::vector< bool > &pairedMinima, const std::vector< bool > &paired1Saddles, const std::vector< bool > &paired2Saddles, const std::vector< bool > &pairedMaxima) const
Print number of pairs, critical cells per dimension & unpaired cells.
void tripletsToPersistencePairs(std::vector< PersistencePair > &pairs, std::vector< bool > &pairedExtrema, std::vector< bool > &pairedSaddles, std::vector< SimplexId > &reps, std::vector< tripletType > &triplets, const SimplexId *const saddlesOrder, const SimplexId *const extremaOrder, const SimplexId pairDim) const
Compute persistence pairs from triplets.
int SimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:22
Persistence pair struct as exported by DiscreteGradient.
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)