TTK
Loading...
Searching...
No Matches
FTMTree_MT_Template.h
Go to the documentation of this file.
1
15
16#pragma once
17
18#include <functional>
19
20#include "FTMTree_MT.h"
21
22#ifdef TTK_ENABLE_OMP_PRIORITY
23#define OPTIONAL_PRIORITY(value) priority(value)
24#else
25#define OPTIONAL_PRIORITY(value)
26#endif
27
28// ----
29// Init
30// ----
31
32namespace ttk {
33 namespace ftm {
34
35 template <class triangulationType>
36 void FTMTree_MT::build(const triangulationType *mesh, const bool ct) {
37 std::string treeString;
38 // Comparator init (template)
39 initComp();
40 switch(mt_data_.treeType) {
41 case TreeType::Join:
42 treeString = "JT";
43 break;
44 case TreeType::Split:
45 treeString = "ST";
46 break;
47 default:
48 treeString = "CT";
49 break;
50 }
51
52 // Build Merge treeString using tasks
53 Timer precomputeTime;
54 int alreadyDone = leafSearch(mesh);
55 printTime(precomputeTime, "leafSearch " + treeString, 3 + alreadyDone);
56
57 Timer buildTime;
58 leafGrowth(mesh);
59 printTime(buildTime, "leafGrowth " + treeString, 3);
60
61 Timer bbTime;
62 trunk(mesh, ct);
63 printTime(bbTime, "trunk " + treeString, 3);
64
65 if(this->getNumberOfNodes() != this->getNumberOfSuperArcs() + 1) {
66 this->printErr(treeString + " not a tree!");
67 }
68
69 // Segmentation
70 if(ct && params_->segm) {
71 Timer segmTime;
73 printTime(segmTime, "segment " + treeString, 3);
74 }
75 }
76
77 // ------------------------------------------------------------------------
78
79 template <class triangulationType>
80 int FTMTree_MT::leafSearch(const triangulationType *mesh) {
81 int ret = 0;
82 // if not already computed by CT
83 if(getNumberOfNodes() == 0) {
84 const auto nbScalars = scalars_->size;
85 const auto chunkSize = getChunkSize();
86 const auto chunkNb = getChunkCount();
87
88 // Extrema extract and launch tasks
89 for(SimplexId chunkId = 0; chunkId < chunkNb; ++chunkId) {
90#ifdef TTK_ENABLE_OPENMP
91#pragma omp task firstprivate(chunkId) OPTIONAL_PRIORITY(isPrior())
92#endif
93 {
94 const SimplexId lowerBound = chunkId * chunkSize;
95 const SimplexId upperBound
96 = std::min(nbScalars, (chunkId + 1) * chunkSize);
97 for(SimplexId v = lowerBound; v < upperBound; ++v) {
98 const auto &neighNumb = mesh->getVertexNeighborNumber(v);
99 valence val = 0;
100
101 for(valence n = 0; n < neighNumb; ++n) {
102 SimplexId neigh{-1};
103 mesh->getVertexNeighbor(v, n, neigh);
104 comp_.vertLower(neigh, v) && ++val;
105 }
106
107 mt_data_.valences[v] = val;
108
109 if(!val) {
110 makeNode(v);
111 }
112 }
113 }
114 }
115
116#ifdef TTK_ENABLE_OPENMP
117#pragma omp taskwait
118#endif
119 } else {
120 ret = 1;
121 }
122
123 // fill leaves
124 const auto &nbLeaves = mt_data_.nodes->size();
125 mt_data_.leaves.resize(nbLeaves);
126 std::iota(mt_data_.leaves.begin(), mt_data_.leaves.end(), 0);
127
128 if(debugLevel_ >= 4) {
129 this->printMsg("found " + std::to_string(nbLeaves) + " leaves");
130 }
131
132 // Reserve Arcs
133 mt_data_.superArcs->reserve(nbLeaves * 2 + 1);
134#ifdef TTK_ENABLE_FTM_TREE_STATS_TIME
135 createVector<ActiveTask>(mt_data_.activeTasksStats);
136 mt_data_.activeTasksStats.resize(nbLeaves * 2 + 1);
137#endif
138
139 return ret;
140 }
141
142 // ------------------------------------------------------------------------
143
144 template <class triangulationType>
145 void FTMTree_MT::leafGrowth(const triangulationType *mesh) {
146 _launchGlobalTime.reStart();
147
148 const auto &nbLeaves = mt_data_.leaves.size();
149
150 // memory allocation here
151 initVectStates(nbLeaves + 2);
152
153 // elevation: backbone only
154 if(nbLeaves == 1) {
155 const SimplexId v = (*mt_data_.nodes)[0].getVertexId();
156 mt_data_.openedNodes[v] = 1;
157 mt_data_.storage.emplace_back(v);
158 mt_data_.ufs[v] = &mt_data_.storage[0];
159 return;
160 }
161
162 mt_data_.activeTasks = nbLeaves;
163 mt_data_.storage.resize(nbLeaves);
164
165 auto comp = [this](const idNode a, const idNode b) {
166#ifdef HIGHER
167 return this->comp_.vertHigher(
168 this->getNode(a)->getVertexId(), this->getNode(b)->getVertexId());
169#else
170 return this->comp_.vertLower(
171 this->getNode(a)->getVertexId(), this->getNode(b)->getVertexId());
172#endif
173 };
174 sort(mt_data_.leaves.begin(), mt_data_.leaves.end(), comp);
175
176 for(idNode n = 0; n < nbLeaves; ++n) {
177 const idNode l = mt_data_.leaves[n];
178 SimplexId v = getNode(l)->getVertexId();
179 // for each node: get vert, create uf and launch
180 mt_data_.storage[n] = AtomicUF{v};
181 mt_data_.ufs[v] = &mt_data_.storage[n];
182
183#ifdef TTK_ENABLE_OPENMP
184#pragma omp task UNTIED() OPTIONAL_PRIORITY(isPrior())
185#endif
186 arcGrowth(mesh, v, n);
187 }
188
189#ifdef TTK_ENABLE_OPENMP
190#pragma omp taskwait
191#endif
192 }
193
194 // ------------------------------------------------------------------------
195
196 template <class triangulationType>
197 void FTMTree_MT::arcGrowth(const triangulationType *mesh,
198 const SimplexId startVert,
199 const SimplexId orig) {
200 // current task id / propag
201
202 // local order (ignore non regular verts)
203 SimplexId localOrder = -1;
204 UF startUF = mt_data_.ufs[startVert]->find();
205 // get or recover states
206 CurrentState *currentState;
207 if(startUF->getNbStates()) {
208 currentState = startUF->getFirstState();
209 } else {
210 const std::size_t currentStateId = mt_data_.states->getNext();
211 currentState = &(*mt_data_.states)[currentStateId];
212 currentState->setStartVert(startVert);
213 startUF->addState(currentState);
214 }
215
216 currentState->addNewVertex(startVert);
217
218 // avoid duplicate processing of startVert
219 bool seenFirst = false;
220
221 // ARC OPENING
222 idNode startNode = getCorrespondingNodeId(startVert);
223 idSuperArc currentArc = openSuperArc(startNode);
224 startUF->addArcToClose(currentArc);
225#ifdef TTK_ENABLE_FTM_TREE_STATS_TIME
226 (*mt_data_.activeTasksStats)[currentArc].begin
227 = _launchGlobalTime.getElapsedTime();
228 (*mt_data_.activeTasksStats)[currentArc].origin = orig;
229#endif
230
231 // TASK PROPAGATION
232 while(!currentState->empty()) {
233 // Next vertex
234
235 SimplexId currentVert = currentState->getNextMinVertex();
236
237 // ignore duplicate
238 if(!isCorrespondingNull(currentVert)
239 && !isCorrespondingNode(currentVert)) {
240 continue;
241 } else {
242 // first node can be duplicate, avoid duplicate process
243 if(currentVert == startVert) {
244 if(!seenFirst) {
245 seenFirst = true;
246 } else {
247 continue;
248 }
249 }
250 }
251
252 // local order to avoid sort
253 mt_data_.visitOrder[currentVert] = localOrder++;
254
255 // Saddle & Last detection + propagation
256 bool isSaddle, isLast;
257 std::tie(isSaddle, isLast) = propagate(mesh, *currentState, startUF);
258
259 // regular propagation
260#ifdef TTK_ENABLE_OPENMP
261#pragma omp atomic write seq_cst
262#endif
263 mt_data_.ufs[currentVert] = startUF;
264
265 // Saddle case
266 if(isSaddle) {
267
268#ifdef TTK_ENABLE_FTM_TREE_STATS_TIME
269 (*mt_data_.activeTasksStats)[currentArc].end
270 = _launchGlobalTime.getElapsedTime();
271#endif
272 // need a node on this vertex
273#ifdef TTK_ENABLE_OPENMP
274#pragma omp atomic write seq_cst
275#endif
276 mt_data_.openedNodes[currentVert] = 1;
277
278 // If last close all and merge
279 if(isLast) {
280 // finish works here
281 closeAndMergeOnSaddle(mesh, currentVert);
282
283 // last task detection
284 idNode remainingTasks;
285#ifdef TTK_ENABLE_OPENMP
286#pragma omp atomic read seq_cst
287#endif
288 remainingTasks = mt_data_.activeTasks;
289 if(remainingTasks == 1) {
290 // only backbone remaining
291 return;
292 }
293
294 // made a node on this vertex
295#ifdef TTK_ENABLE_OPENMP
296#pragma omp atomic write seq_cst
297#endif
298 mt_data_.openedNodes[currentVert] = 0;
299
300 // recursively continue
301#ifdef TTK_ENABLE_OPENMP
302#pragma omp taskyield
303#endif
304 arcGrowth(mesh, currentVert, orig);
305 } else {
306 // Active tasks / threads
307#ifdef TTK_ENABLE_OPENMP
308#pragma omp atomic update seq_cst
309#endif
311 }
312
313 // stop at saddle
314 return;
315 }
316
317 if(currentVert != startVert) {
318 updateCorrespondingArc(currentVert, currentArc);
319 }
320 getSuperArc(currentArc)->setLastVisited(currentVert);
321
322 } // end wile propagation
323
324 // close root
325 const SimplexId closeVert = getSuperArc(currentArc)->getLastVisited();
326 bool existCloseNode = isCorrespondingNode(closeVert);
327 idNode closeNode = (existCloseNode) ? getCorrespondingNodeId(closeVert)
328 : makeNode(closeVert);
329 closeSuperArc(currentArc, closeNode);
330 getSuperArc(currentArc)->decrNbSeen();
331 idNode rootPos = mt_data_.roots->getNext();
332 (*mt_data_.roots)[rootPos] = closeNode;
333
334#ifdef TTK_ENABLE_FTM_TREE_STATS_TIME
335 mt_data_.activeTasksStats[currentArc].end
336 = _launchGlobalTime.getElapsedTime();
337#endif
338 }
339
340 // ------------------------------------------------------------------------
341
342 template <class triangulationType>
343 std::tuple<bool, bool> FTMTree_MT::propagate(const triangulationType *mesh,
344 CurrentState &currentState,
345 UF curUF) {
346 bool becameSaddle = false, isLast = false;
347 const auto nbNeigh = mesh->getVertexNeighborNumber(currentState.vertex);
348 valence decr = 0;
349
350 // once for all
351 auto *curUFF = curUF->find();
352
353 // propagation / is saddle
354 for(valence n = 0; n < nbNeigh; ++n) {
355 SimplexId neigh{-1};
356 mesh->getVertexNeighbor(currentState.vertex, n, neigh);
357
358 if(comp_.vertLower(neigh, currentState.vertex)) {
359 UF neighUF = mt_data_.ufs[neigh];
360
361 // is saddle
362 if(!neighUF || neighUF->find() != curUFF) {
363 becameSaddle = true;
364 } else if(neighUF) {
365 ++decr;
366 }
367
368 } else {
369 if(!mt_data_.propagation[neigh]
370 || mt_data_.propagation[neigh]->find() != curUFF) {
371 currentState.addNewVertex(neigh);
372 mt_data_.propagation[neigh] = curUFF;
373 }
374 }
375 }
376
377 // is last
378 valence oldVal;
379 valence &tmp = mt_data_.valences[currentState.vertex];
380#ifdef TTK_ENABLE_OPENMP
381#pragma omp atomic capture
382#endif
383 {
384 oldVal = tmp;
385 tmp -= decr;
386 }
387 if(oldVal == decr) {
388 isLast = true;
389 }
390
391 return std::make_tuple(becameSaddle, isLast);
392 }
393
394 // ------------------------------------------------------------------------
395
396 template <class triangulationType>
397 SimplexId FTMTree_MT::trunk(const triangulationType *mesh, const bool ct) {
398 Timer bbTimer;
399
400 std::vector<SimplexId> trunkVerts;
401 const auto &nbScalars = scalars_->size;
402
403 // trunkVerts
404 trunkVerts.reserve(std::max(SimplexId{10}, nbScalars / 500));
405 for(SimplexId v = 0; v < nbScalars; ++v) {
406 if(mt_data_.openedNodes[v]) {
407 if(this->isCorrespondingNode(v)) {
408 // parallel leafGrowth can partially build the trunk,
409 // filter out the saddles with an upward arc
410 const auto node{this->getNode(this->getCorrespondingNodeId(v))};
411 if(node->getNumberOfUpSuperArcs() == 0) {
412 trunkVerts.emplace_back(v);
413 }
414 } else {
415 trunkVerts.emplace_back(v);
416 }
417 }
418 }
419 sort(trunkVerts.begin(), trunkVerts.end(), comp_.vertLower);
420 for(const SimplexId v : trunkVerts) {
421 closeOnBackBone(mesh, v);
422 }
423
424 // Arcs
425 const auto &nbNodes = trunkVerts.size();
426 for(idNode n = 1; n < nbNodes; ++n) {
427 idSuperArc na = makeSuperArc(getCorrespondingNodeId(trunkVerts[n - 1]),
428 getCorrespondingNodeId(trunkVerts[n]));
429 getSuperArc(na)->setLastVisited(trunkVerts[n]);
430 }
431
432 if(!nbNodes) {
433 return 0;
434 }
435 const idSuperArc lastArc
436 = openSuperArc(getCorrespondingNodeId(trunkVerts[nbNodes - 1]));
437
438 // Root (close last arc)
439 // if several CC still the backbone is only in one.
440 // But the root may not be the max node of the whole dataset: TODO
441 const idNode rootNode
442 = makeNode(scalars_->sortedVertices[(isJT()) ? scalars_->size - 1 : 0]);
443 closeSuperArc(lastArc, rootNode);
444 getSuperArc(lastArc)->setLastVisited(getNode(rootNode)->getVertexId());
445
446 printTime(bbTimer, "trunk seq.", 4);
447 bbTimer.reStart();
448
449 // Segmentation
450 SimplexId begin, stop, processed;
451 std::tie(begin, stop) = getBoundsFromVerts(trunkVerts);
452 if(ct) {
453 processed = trunkCTSegmentation(trunkVerts, begin, stop);
454 } else {
455 processed = trunkSegmentation(trunkVerts, begin, stop);
456 }
457 printTime(bbTimer, "trunk para.", 4);
458
459 return processed;
460 }
461
462 // ------------------------------------------------------------------------
463
464 template <class triangulationType>
465 void FTMTree_MT::closeAndMergeOnSaddle(const triangulationType *mesh,
466 SimplexId saddleVert) {
467 idNode closeNode = makeNode(saddleVert);
468
469 // Union of the UF coming here (merge propagation and closing arcs)
470 const auto &nbNeigh = mesh->getVertexNeighborNumber(saddleVert);
471 for(valence n = 0; n < nbNeigh; ++n) {
472 SimplexId neigh{-1};
473 mesh->getVertexNeighbor(saddleVert, n, neigh);
474
475 if(comp_.vertLower(neigh, saddleVert)) {
476 if(mt_data_.ufs[neigh]->find() != mt_data_.ufs[saddleVert]->find()) {
477 mt_data_.ufs[saddleVert] = AtomicUF::makeUnion(
478 mt_data_.ufs[saddleVert], mt_data_.ufs[neigh]);
479 }
480 }
481 }
482
483 // close arcs on this node
484 closeArcsUF(closeNode, mt_data_.ufs[saddleVert]);
485
486 mt_data_.ufs[saddleVert]->find()->mergeStates();
487 mt_data_.ufs[saddleVert]->find()->setExtrema(saddleVert);
488 }
489
490 // ------------------------------------------------------------------------
491
492 template <class triangulationType>
493 void FTMTree_MT::closeOnBackBone(const triangulationType *mesh,
494 SimplexId saddleVert) {
495 idNode closeNode = makeNode(saddleVert);
496
497 // Union of the UF coming here (merge propagation and closing arcs)
498 const auto &nbNeigh = mesh->getVertexNeighborNumber(saddleVert);
499 for(valence n = 0; n < nbNeigh; ++n) {
500 SimplexId neigh{-1};
501 mesh->getVertexNeighbor(saddleVert, n, neigh);
502
503 if(comp_.vertLower(neigh, saddleVert)) {
504 if(mt_data_.ufs[neigh]
505 && mt_data_.ufs[neigh]->find()
506 != mt_data_.ufs[saddleVert]->find()) {
507 mt_data_.ufs[saddleVert] = AtomicUF::makeUnion(
508 mt_data_.ufs[saddleVert], mt_data_.ufs[neigh]);
509 }
510 }
511 }
512
513 // close arcs on this node
514 closeArcsUF(closeNode, mt_data_.ufs[saddleVert]);
515 }
516
517 // ------------------------------------------------------------------------
518
519 template <typename scalarType>
521
522 const auto nbVertices = scalars_->size;
523 scalars_->sortedVertices.resize(nbVertices);
524
525#ifdef TTK_ENABLE_OPENMP
526#pragma omp parallel for
527#endif
528 for(SimplexId i = 0; i < nbVertices; i++) {
529 scalars_->sortedVertices[scalars_->offsets[i]] = i;
530 }
531 }
532
533 } // namespace ftm
534} // namespace ttk
535// Process
int debugLevel_
Definition: Debug.h:379
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
int printErr(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition: Debug.h:149
void reStart()
Definition: Timer.h:21
void addArcToClose(idSuperArc a)
Definition: FTMAtomicUF.h:106
CurrentState * getFirstState()
Definition: FTMAtomicUF.h:58
void addState(CurrentState *s)
Definition: FTMAtomicUF.h:102
size_t getNbStates() const
Definition: FTMAtomicUF.h:82
static AtomicUF * makeUnion(AtomicUF *uf0, AtomicUF *uf1)
Definition: FTMAtomicUF.h:142
AtomicUF * find()
Definition: FTMAtomicUF.h:38
bool isCorrespondingNode(const SimplexId val) const
Definition: FTMTree_MT.h:449
std::shared_ptr< Params > params_
Definition: FTMTree_MT.h:86
void closeOnBackBone(const triangulationType *mesh, SimplexId saddleVert)
SuperArc * getSuperArc(idSuperArc i)
Definition: FTMTree_MT.h:367
idNode getNumberOfNodes() const
Definition: FTMTree_MT.h:389
SimplexId getChunkSize(const SimplexId nbVerts=-1, const SimplexId nbtasks=100) const
Definition: FTMTree_MT.h:819
bool isCorrespondingNull(const SimplexId val) const
Definition: FTMTree_MT.h:453
void sortInput()
if sortedVertices_ is null, define and fill it
std::tuple< SimplexId, SimplexId > getBoundsFromVerts(const std::vector< SimplexId > &nodes) const
Definition: FTMTree_MT.cpp:353
int leafSearch(const triangulationType *mesh)
void closeArcsUF(idNode closeNode, UF uf)
Definition: FTMTree_MT.cpp:264
void closeAndMergeOnSaddle(const triangulationType *mesh, SimplexId saddleVert)
virtual SimplexId trunkSegmentation(const std::vector< SimplexId > &pendingNodesVerts, const SimplexId begin, const SimplexId stop)
std::tuple< bool, bool > propagate(const triangulationType *mesh, CurrentState &currentState, UF curUF)
idNode getCorrespondingNodeId(const SimplexId val) const
Definition: FTMTree_MT.h:459
idSuperArc getNumberOfSuperArcs() const
Definition: FTMTree_MT.h:363
void build(const triangulationType *mesh, const bool ct)
Compute the merge.
int printTime(Timer &t, const std::string &s, const int debugLevel=2) const
Definition: FTMTree_MT.cpp:727
SimplexId getChunkCount(const SimplexId nbVerts=-1, const SimplexId nbTasks=100) const
Definition: FTMTree_MT.h:832
SimplexId trunkCTSegmentation(const std::vector< SimplexId > &pendingNodesVerts, const SimplexId begin, const SimplexId stop)
Definition: FTMTree_MT.cpp:931
idSuperArc openSuperArc(idNode downNodeId)
Definition: FTMTree_MT.cpp:627
SimplexId trunk(const triangulationType *mesh, const bool ct)
idNode makeNode(SimplexId vertexId, SimplexId linked=nullVertex)
Definition: FTMTree_MT.cpp:473
void buildSegmentation()
use vert2tree to compute the segmentation of the fresh built merge tree.
Definition: FTMTree_MT.cpp:102
void leafGrowth(const triangulationType *mesh)
idSuperArc makeSuperArc(idNode downNodeId, idNode upNodeId)
Definition: FTMTree_MT.cpp:499
bool isJT() const
Definition: FTMTree_MT.h:289
Node * getNode(idNode nodeId)
Definition: FTMTree_MT.h:393
void closeSuperArc(idSuperArc superArcId, idNode upNodeId)
Definition: FTMTree_MT.cpp:271
std::shared_ptr< Scalars > scalars_
Definition: FTMTree_MT.h:87
void initVectStates(const SimplexId nbLeaves)
Definition: FTMTree_MT.h:202
Comparison comp_
Definition: FTMTree_MT.h:91
void updateCorrespondingArc(const SimplexId vert, const idSuperArc arc)
Definition: FTMTree_MT.h:493
void arcGrowth(const triangulationType *mesh, const SimplexId startVert, const SimplexId orig)
SimplexId getVertexId() const
Definition: FTMNode.h:54
SimplexId getLastVisited() const
Definition: FTMSuperArc.h:86
void setLastVisited(SimplexId vertId)
Definition: FTMSuperArc.h:90
long unsigned int idSuperArc
SuperArc index in vect_superArcs_.
Definition: FTMDataTypes.h:37
SimplexId valence
for vertex up/down valence
Definition: FTMDataTypes.h:63
unsigned int idNode
Node index in vect_nodes_.
Definition: FTMDataTypes.h:39
The Topology ToolKit.
int SimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:22
void setStartVert(const SimplexId v)
Definition: FTMStructures.h:81
SimplexId getNextMinVertex()
Definition: FTMStructures.h:85
void addNewVertex(const SimplexId v)
Definition: FTMStructures.h:91
std::shared_ptr< FTMAtomicVector< SuperArc > > superArcs
Definition: FTMTree_MT.h:46
std::vector< AtomicUF > storage
Definition: FTMTree_MT.h:57
std::shared_ptr< FTMAtomicVector< CurrentState > > states
Definition: FTMTree_MT.h:60
std::shared_ptr< FTMAtomicVector< Node > > nodes
Definition: FTMTree_MT.h:47
std::vector< SimplexId > visitOrder
Definition: FTMTree_MT.h:53
std::vector< UF > propagation
Definition: FTMTree_MT.h:59
std::vector< UF > ufs
Definition: FTMTree_MT.h:58
std::shared_ptr< FTMAtomicVector< idNode > > roots
Definition: FTMTree_MT.h:48
std::vector< valence > valences
Definition: FTMTree_MT.h:62
TreeType treeType
Definition: FTMTree_MT.h:43
std::vector< idNode > leaves
Definition: FTMTree_MT.h:49
std::vector< char > openedNodes
Definition: FTMTree_MT.h:64