TTK
Loading...
Searching...
No Matches
PersistenceDiagramBarycenter.cpp
Go to the documentation of this file.
2
4 std::vector<DiagramType> &intermediateDiagrams,
5 DiagramType &barycenter,
6 std::vector<std::vector<std::vector<MatchingType>>> &all_matchings) {
7
8 Timer tm;
9
10 printMsg("Computing Barycenter of " + std::to_string(numberOfInputs_)
11 + " diagrams.");
12
13 std::vector<DiagramType> data_min(numberOfInputs_);
14 std::vector<DiagramType> data_sad(numberOfInputs_);
15 std::vector<DiagramType> data_max(numberOfInputs_);
16
17 std::vector<std::vector<int>> data_min_idx(numberOfInputs_);
18 std::vector<std::vector<int>> data_sad_idx(numberOfInputs_);
19 std::vector<std::vector<int>> data_max_idx(numberOfInputs_);
20
21 bool do_min = false;
22 bool do_sad = false;
23 bool do_max = false;
24
25 // Create diagrams for min, saddle and max persistence pairs
26 for(int i = 0; i < numberOfInputs_; i++) {
27 DiagramType &CTDiagram = intermediateDiagrams[i];
28
29 for(size_t j = 0; j < CTDiagram.size(); ++j) {
30 PersistencePair const &t = CTDiagram[j];
31
32 ttk::CriticalType const nt1 = t.birth.type;
33 ttk::CriticalType const nt2 = t.death.type;
34
35 double const dt = t.persistence();
36 // if (abs<double>(dt) < zeroThresh) continue;
37 if(dt > 0) {
40 data_max[i].push_back(t);
41 data_max_idx[i].push_back(j);
42 do_max = true;
43 } else {
46 data_max[i].push_back(t);
47 data_max_idx[i].push_back(j);
48 do_max = true;
49 }
52 data_min[i].push_back(t);
53 data_min_idx[i].push_back(j);
54 do_min = true;
55 }
59 && nt2 == ttk::CriticalType::Saddle1)) {
60 data_sad[i].push_back(t);
61 data_sad_idx[i].push_back(j);
62 do_sad = true;
63 }
64 }
65 }
66 }
67 }
68
69 DiagramType barycenter_min;
70 DiagramType barycenter_sad;
71 DiagramType barycenter_max;
72
73 std::vector<std::vector<MatchingType>> matching_min, matching_sad,
74 matching_max;
75
76 double total_cost = 0, min_cost = 0, sad_cost = 0, max_cost = 0;
77 /*omp_set_num_threads(1);
78 #ifdef TTK_ENABLE_OPENMP
79 #pragma omp parallel sections
80 #endif
81 {
82 #ifdef TTK_ENABLE_OPENMP
83 #pragma omp section
84 #endif
85 {*/
86 if(do_min) {
87 printMsg("Computing Minima barycenter...");
88 PDBarycenter bary_min{};
90 bary_min.setWasserstein(wasserstein_);
91 bary_min.setNumberOfInputs(numberOfInputs_);
92 bary_min.setDiagramType(0);
93 bary_min.setUseProgressive(use_progressive_);
94 bary_min.setGeometricalFactor(alpha_);
95 bary_min.setDebugLevel(debugLevel_);
96 bary_min.setDeterministic(deterministic_);
97 bary_min.setLambda(lambda_);
98 bary_min.setMethod(method_);
99 bary_min.setEarlyStoppage(early_stoppage_);
100 bary_min.setEpsilonDecreases(epsilon_decreases_);
101 bary_min.setReinitPrices(reinit_prices_);
102 bary_min.setNonMatchingWeight(nonMatchingWeight_);
103 bary_min.setDiagrams(&data_min);
104 bary_min.setDeltaLim(delta_lim_);
105 matching_min = bary_min.execute(barycenter_min);
106 min_cost = bary_min.getCost();
107 total_cost += min_cost;
108 }
109 /*}
110
111 #ifdef TTK_ENABLE_OPENMP
112 #pragma omp section
113 #endif
114 {*/
115 if(do_sad) {
116 printMsg("Computing Saddles barycenter...");
117 PDBarycenter bary_sad{};
119 bary_sad.setWasserstein(wasserstein_);
120 bary_sad.setNumberOfInputs(numberOfInputs_);
121 bary_sad.setDiagramType(1);
122 bary_sad.setUseProgressive(use_progressive_);
123 bary_sad.setGeometricalFactor(alpha_);
124 bary_sad.setLambda(lambda_);
125 bary_sad.setDebugLevel(debugLevel_);
126 bary_sad.setMethod(method_);
127 bary_sad.setEarlyStoppage(early_stoppage_);
128 bary_sad.setEpsilonDecreases(epsilon_decreases_);
129 bary_sad.setDeterministic(deterministic_);
130 bary_sad.setReinitPrices(reinit_prices_);
131 bary_sad.setNonMatchingWeight(nonMatchingWeight_);
132 bary_sad.setDiagrams(&data_sad);
133 bary_sad.setDeltaLim(delta_lim_);
134 matching_sad = bary_sad.execute(barycenter_sad);
135 sad_cost = bary_sad.getCost();
136 total_cost += sad_cost;
137 }
138 /*}
139
140 #ifdef TTK_ENABLE_OPENMP
141 #pragma omp section
142 #endif
143 {*/
144 if(do_max) {
145 printMsg("Computing Maxima barycenter...");
146 PDBarycenter bary_max{};
148 bary_max.setWasserstein(wasserstein_);
149 bary_max.setNumberOfInputs(numberOfInputs_);
150 bary_max.setDiagramType(2);
151 bary_max.setUseProgressive(use_progressive_);
152 bary_max.setGeometricalFactor(alpha_);
153 bary_max.setLambda(lambda_);
154 bary_max.setMethod(method_);
155 bary_max.setDebugLevel(debugLevel_);
156 bary_max.setEarlyStoppage(early_stoppage_);
157 bary_max.setDeterministic(deterministic_);
158 bary_max.setEpsilonDecreases(epsilon_decreases_);
159 bary_max.setReinitPrices(reinit_prices_);
160 bary_max.setNonMatchingWeight(nonMatchingWeight_);
161 bary_max.setDiagrams(&data_max);
162 bary_max.setDeltaLim(delta_lim_);
163 matching_max = bary_max.execute(barycenter_max);
164 max_cost = bary_max.getCost();
165 total_cost += max_cost;
166 }
167 //}
168 //}
169
170 // Reconstruct matchings
171 all_matchings.resize(1);
172 all_matchings[0].resize(numberOfInputs_);
173 for(int i = 0; i < numberOfInputs_; i++) {
174
175 if(do_min) {
176 for(size_t j = 0; j < matching_min[i].size(); j++) {
177 MatchingType t = matching_min[i][j];
178 int const bidder_id = std::get<0>(t);
179 std::get<0>(t) = data_min_idx[i][bidder_id];
180 if(std::get<1>(t) < 0) {
181 std::get<1>(t) = -1;
182 }
183 all_matchings[0][i].push_back(t);
184 }
185 }
186
187 if(do_sad) {
188 for(size_t j = 0; j < matching_sad[i].size(); j++) {
189 MatchingType t = matching_sad[i][j];
190 int const bidder_id = std::get<0>(t);
191 std::get<0>(t) = data_sad_idx[i][bidder_id];
192 if(std::get<1>(t) >= 0) {
193 std::get<1>(t) = std::get<1>(t) + barycenter_min.size();
194 } else {
195 std::get<1>(t) = -1;
196 }
197 all_matchings[0][i].push_back(t);
198 }
199 }
200
201 if(do_max) {
202 for(size_t j = 0; j < matching_max[i].size(); j++) {
203 MatchingType t = matching_max[i][j];
204 int const bidder_id = std::get<0>(t);
205 std::get<0>(t) = data_max_idx[i][bidder_id];
206 if(std::get<1>(t) >= 0) {
207 std::get<1>(t)
208 = std::get<1>(t) + barycenter_min.size() + barycenter_sad.size();
209 } else {
210 std::get<1>(t) = -1;
211 }
212 all_matchings[0][i].push_back(t);
213 }
214 }
215 }
216 // Reconstruct barcenter
217 for(size_t j = 0; j < barycenter_min.size(); j++) {
218 const auto &dt = barycenter_min[j];
219 barycenter.push_back(dt);
220 }
221 for(size_t j = 0; j < barycenter_sad.size(); j++) {
222 const auto &dt = barycenter_sad[j];
223 barycenter.push_back(dt);
224 }
225 for(size_t j = 0; j < barycenter_max.size(); j++) {
226 const auto &dt = barycenter_max[j];
227 barycenter.push_back(dt);
228 }
229
230 // Recreate 3D critical coordinates of barycentric points
231 std::vector<int> number_of_matchings_for_point(barycenter.size());
232 std::vector<float> cords_x1(barycenter.size());
233 std::vector<float> cords_y1(barycenter.size());
234 std::vector<float> cords_z1(barycenter.size());
235 std::vector<float> cords_x2(barycenter.size());
236 std::vector<float> cords_y2(barycenter.size());
237 std::vector<float> cords_z2(barycenter.size());
238 for(unsigned i = 0; i < barycenter.size(); i++) {
239 number_of_matchings_for_point[i] = 0;
240 cords_x1[i] = 0;
241 cords_y1[i] = 0;
242 cords_z1[i] = 0;
243 cords_x2[i] = 0;
244 cords_y2[i] = 0;
245 cords_z2[i] = 0;
246 }
247
248 for(unsigned i = 0; i < all_matchings[0].size(); i++) {
249 DiagramType &CTDiagram = intermediateDiagrams[i];
250 for(unsigned j = 0; j < all_matchings[0][i].size(); j++) {
251 MatchingType t = all_matchings[0][i][j];
252 int const bidder_id = std::get<0>(t);
253 int const bary_id = std::get<1>(t);
254
255 const auto &bidder = CTDiagram[bidder_id];
256 number_of_matchings_for_point[bary_id] += 1;
257 cords_x1[bary_id] += bidder.birth.coords[0];
258 cords_y1[bary_id] += bidder.birth.coords[1];
259 cords_z1[bary_id] += bidder.birth.coords[2];
260 cords_x2[bary_id] += bidder.death.coords[0];
261 cords_y2[bary_id] += bidder.death.coords[1];
262 cords_z2[bary_id] += bidder.death.coords[2];
263 }
264 }
265
266 for(unsigned i = 0; i < barycenter.size(); i++) {
267 if(number_of_matchings_for_point[i] > 0) {
268 barycenter[i].birth.coords[0]
269 = cords_x1[i] / number_of_matchings_for_point[i];
270 barycenter[i].birth.coords[1]
271 = cords_y1[i] / number_of_matchings_for_point[i];
272 barycenter[i].birth.coords[2]
273 = cords_z1[i] / number_of_matchings_for_point[i];
274 barycenter[i].death.coords[0]
275 = cords_x2[i] / number_of_matchings_for_point[i];
276 barycenter[i].death.coords[1]
277 = cords_y2[i] / number_of_matchings_for_point[i];
278 barycenter[i].death.coords[2]
279 = cords_z2[i] / number_of_matchings_for_point[i];
280 }
281 }
282
283 printMsg("Min-saddle cost : " + std::to_string(min_cost));
284 printMsg("Saddle-saddle cost : " + std::to_string(sad_cost));
285 printMsg("Saddle-max cost : " + std::to_string(max_cost));
286 printMsg("Total cost : " + std::to_string(total_cost));
287 printMsg("Complete", 1, tm.getElapsedTime(), threadNumber_);
288}
virtual int setThreadNumber(const int threadNumber)
Definition BaseClass.h:80
int debugLevel_
Definition Debug.h:379
void execute(std::vector< DiagramType > &intermediateDiagrams, DiagramType &barycenter, std::vector< std::vector< std::vector< MatchingType > > > &all_matchings)
double getElapsedTime()
Definition Timer.h:15
std::string to_string(__int128)
Definition ripserpy.cpp:99
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.
double persistence() const
Return the topological persistence of the pair.
printMsg(debug::output::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/|__ _|"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)