TTK
Loading...
Searching...
No Matches
PersistenceDiagram.h
Go to the documentation of this file.
1
159
160#pragma once
161
162// base code includes
163#include <ApproximateTopology.h>
166#include <FTMTreePP.h>
169#include <ProgressiveTopology.h>
170#include <Triangulation.h>
171#include <psort.h>
172
173namespace ttk {
174
175 namespace persistenceSort {
176 inline bool comp(const PersistencePair a, const PersistencePair b) {
177 return a.birth.offset < b.birth.offset;
178 };
179
180 inline bool oppositeComp(const PersistencePair a, const PersistencePair b) {
181 return a.birth.offset > b.birth.offset;
182 };
183 } // namespace persistenceSort
184
189 class PersistenceDiagram : virtual public Debug {
190
191 public:
200
202
203 inline void setBackend(const BACKEND be) {
204 this->BackEnd = be;
205 }
206
207 inline void setComputeMinSad(const bool data) {
208 this->dms_.setComputeMinSad(data);
209#if defined(TTK_ENABLE_MPI) && defined(TTK_ENABLE_OPENMP)
210 this->dmsMPI_.setComputeMinSad(data);
211#endif
212 }
213 inline void setComputeSadSad(const bool data) {
214 this->dms_.setComputeSadSad(data);
215#if defined(TTK_ENABLE_MPI) && defined(TTK_ENABLE_OPENMP)
216 this->dmsMPI_.setComputeSadSad(data);
217#endif
218 }
219 inline void setComputeSadMax(const bool data) {
220 this->dms_.setComputeSadMax(data);
221#if defined(TTK_ENABLE_MPI) && defined(TTK_ENABLE_OPENMP)
222 this->dmsMPI_.setComputeSadMax(data);
223#endif
224 }
225
230 template <typename scalarType, typename triangulationType>
231 void
232 augmentPersistenceDiagram(std::vector<PersistencePair> &persistencePairs,
233 const scalarType *const scalars,
234 const triangulationType *triangulation);
235
237 ftm::TreeType treeType,
238 const SimplexId vertexId) const;
239
240 void sortPersistenceDiagram(std::vector<PersistencePair> &diagram,
241 const SimplexId *const offsets) const;
242
243 template <typename scalarType>
245 ftm::FTMTreePP &tree,
246 const std::vector<
247 std::tuple<ttk::SimplexId, ttk::SimplexId, scalarType, bool>> &pairs,
248 std::vector<PersistencePair> &diagram) const;
249
258 template <typename scalarType, class triangulationType>
259 int execute(std::vector<PersistencePair> &CTDiagram,
260 const scalarType *inputScalars,
261 const size_t scalarsMTime,
262 const SimplexId *inputOffsets,
263 const triangulationType *triangulation,
264 const std::vector<bool> *updateMask = nullptr);
265
266 template <typename scalarType, class triangulationType>
267 int executeFTM(std::vector<PersistencePair> &CTDiagram,
268 const scalarType *inputScalars,
269 const SimplexId *inputOffsets,
270 const triangulationType *triangulation);
271
272 template <class triangulationType>
273 int executeProgressiveTopology(std::vector<PersistencePair> &CTDiagram,
274 const SimplexId *inputOffsets,
275 const triangulationType *triangulation);
276 template <typename scalarType, class triangulationType>
277 int executeApproximateTopology(std::vector<PersistencePair> &CTDiagram,
278 const scalarType *inputScalars,
279 const triangulationType *triangulation);
280
281 template <class triangulationType>
282 int executePersistentSimplex(std::vector<PersistencePair> &CTDiagram,
283 const SimplexId *inputOffsets,
284 const triangulationType *triangulation);
285
286#ifdef TTK_ENABLE_MPI
287 template <typename scalarType, class triangulationType>
288 int executeDiscreteMorseSandwichMPI(std::vector<PersistencePair> &CTDiagram,
289 const scalarType *inputScalars,
290 const size_t scalarsMTime,
291 const SimplexId *inputOffsets,
292 const triangulationType *triangulation);
293#endif
294
295 template <typename scalarType, class triangulationType>
296 int executeDiscreteMorseSandwich(std::vector<PersistencePair> &CTDiagram,
297 const scalarType *inputScalars,
298 const size_t scalarsMTime,
299 const SimplexId *inputOffsets,
300 const triangulationType *triangulation,
301 const std::vector<bool> *updateMask
302 = nullptr);
303
304 template <class triangulationType>
305 void checkProgressivityRequirement(const triangulationType *triangulation);
306
307 template <class triangulationType>
308 void checkManifold(const triangulationType *const triangulation);
309
310 inline void
312 if(triangulation) {
313 triangulation->preconditionBoundaryVertices();
314 if(this->BackEnd == BACKEND::FTM
317 contourTree_.setDebugLevel(debugLevel_);
318 contourTree_.setThreadNumber(threadNumber_);
319 contourTree_.preconditionTriangulation(triangulation);
320 }
322 dms_.setDebugLevel(debugLevel_);
323 dms_.setThreadNumber(threadNumber_);
324 dms_.preconditionTriangulation(triangulation);
325 triangulation->preconditionManifold();
326 }
327#if defined(TTK_ENABLE_MPI) && defined(TTK_ENABLE_OPENMP)
329 dmsMPI_.setDebugLevel(debugLevel_);
330 dmsMPI_.setThreadNumber(threadNumber_);
331 dmsMPI_.preconditionTriangulation(triangulation);
332 triangulation->preconditionManifold();
333 }
334#endif
338 psp_.preconditionTriangulation(triangulation);
339 }
340 }
341 }
342
343 inline void setOutputMonotonyOffsets(void *data) {
345 }
346 inline void setOutputOffsets(void *data) {
347 outputOffsets_ = data;
348 }
349 inline void setOutputScalars(void *data) {
350 outputScalars_ = data;
351 }
352 inline void setDeltaApproximate(double data) {
353 approxT_.setDelta(data);
354 }
355
356 protected:
357 bool IgnoreBoundary{false};
362#if defined(TTK_ENABLE_MPI) && defined(TTK_ENABLE_OPENMP)
363 DiscreteMorseSandwichMPI dmsMPI_{};
364#endif
365 // int BackEnd{0};
367 // progressivity
370
373 bool UseTasks{true};
374 bool IsResumable{false};
375 double TimeLimit{};
376
377 // Approximation
381 double Epsilon;
382 };
383} // namespace ttk
384
385template <typename scalarType>
387 ftm::FTMTreePP &tree,
388 const std::vector<
389 std::tuple<ttk::SimplexId, ttk::SimplexId, scalarType, bool>> &pairs,
390 std::vector<PersistencePair> &diagram) const {
391
392 const ttk::SimplexId numberOfPairs = pairs.size();
393 diagram.resize(numberOfPairs);
394 for(ttk::SimplexId i = 0; i < numberOfPairs; ++i) {
395 const ttk::SimplexId v0 = std::get<0>(pairs[i]);
396 const ttk::SimplexId v1 = std::get<1>(pairs[i]);
397 const bool type = std::get<3>(pairs[i]);
398
399 if(type == true) {
400 diagram[i] = PersistencePair{
402 v0,
403 {},
404 {},
405 {},
408 v1,
409 {},
410 {},
411 {},
413 0, true};
414 } else {
415 diagram[i] = PersistencePair{
417 v1,
418 {},
419 {},
420 {},
423 v0,
424 {},
425 {},
426 {},
428 2, true};
429 }
430 }
431
432 diagram.back().isFinite = false; // global extrema pair is infinite
433
434 return 0;
435}
436
437template <typename scalarType, typename triangulationType>
439 std::vector<PersistencePair> &persistencePairs,
440 const scalarType *const scalars,
441 const triangulationType *triangulation) {
442
443#ifdef TTK_ENABLE_OPENMP
444#pragma omp parallel for num_threads(threadNumber_)
445#endif // TTK_ENABLE_OPENMP
446 for(std::size_t i = 0; i < persistencePairs.size(); ++i) {
447 auto &pair{persistencePairs[i]};
448 triangulation->getVertexPoint(pair.birth.id, pair.birth.coords[0],
449 pair.birth.coords[1], pair.birth.coords[2]);
450 pair.birth.sfValue = scalars[pair.birth.id];
451 triangulation->getVertexPoint(pair.death.id, pair.death.coords[0],
452 pair.death.coords[1], pair.death.coords[2]);
453 pair.death.sfValue = scalars[pair.death.id];
454 }
455}
456
457template <typename scalarType, class triangulationType>
458int ttk::PersistenceDiagram::execute(std::vector<PersistencePair> &CTDiagram,
459 const scalarType *inputScalars,
460 const size_t scalarsMTime,
461 const SimplexId *inputOffsets,
462 const triangulationType *triangulation,
463 const std::vector<bool> *updateMask) {
464
466
467 checkProgressivityRequirement(triangulation);
468 checkManifold(triangulation);
469
470 Timer tm{};
471
472 switch(BackEnd) {
474 executePersistentSimplex(CTDiagram, inputOffsets, triangulation);
475 break;
477 executeDiscreteMorseSandwich(CTDiagram, inputScalars, scalarsMTime,
478 inputOffsets, triangulation, updateMask);
479 break;
481 executeProgressiveTopology(CTDiagram, inputOffsets, triangulation);
482 break;
484 executeApproximateTopology(CTDiagram, inputScalars, triangulation);
485 break;
486 case BACKEND::FTM:
487 executeFTM(CTDiagram, inputScalars, inputOffsets, triangulation);
488 break;
490#ifdef TTK_ENABLE_MPI
491 executeDiscreteMorseSandwichMPI(
492 CTDiagram, inputScalars, scalarsMTime, inputOffsets, triangulation);
493#else
494 this->printWrn("TTK is not compiled with MPI. Running sequentially.");
495 this->printWrn("If you want to run TTK with MPI, compile it with "
496 "TTK_ENABLE_MPI to ON.");
498 CTDiagram, inputScalars, scalarsMTime, inputOffsets, triangulation);
499#endif
500 break;
501 default:
502 printErr("No method was selected");
503 }
504
505 this->printMsg("Complete", 1.0, tm.getElapsedTime(), this->threadNumber_);
506#ifdef TTK_ENABLE_MPI
507 if(!isRunningWithMPI()) {
508#endif
509 // augment persistence pairs with meta-data
510 augmentPersistenceDiagram(CTDiagram, inputScalars, triangulation);
511
512 // finally sort the diagram
513 sortPersistenceDiagram(CTDiagram, inputOffsets);
514#ifdef TTK_ENABLE_MPI
515 }
516#endif
517
519
520 return 0;
521}
522
523template <class triangulationType>
525 std::vector<PersistencePair> &CTDiagram,
526 const SimplexId *inputOffsets,
527 const triangulationType *triangulation) {
528
529 Timer const tm{};
530 const auto dim = triangulation->getDimensionality();
531
532 std::vector<ttk::PersistentSimplexPairs::PersistencePair> pairs{};
533
534 psp_.setDebugLevel(this->debugLevel_);
535 psp_.setThreadNumber(this->threadNumber_);
536 psp_.computePersistencePairs(pairs, inputOffsets, *triangulation);
537 dms_.setInputOffsets(inputOffsets);
538
539 // convert PersistentSimplex pairs (with critical cells id) to PL
540 // pairs (with vertices id)
541
542#ifdef TTK_ENABLE_OPENMP
543#pragma omp parallel for num_threads(threadNumber_)
544#endif // TTK_ENABLE_OPENMP
545 for(size_t i = 0; i < pairs.size(); ++i) {
546 auto &pair{pairs[i]};
547 if(pair.type > 0) {
548 pair.birth = dms_.getCellGreaterVertex(
549 Cell{pair.type, pair.birth}, *triangulation);
550 }
551 if(pair.death != -1) {
552 pair.death = dms_.getCellGreaterVertex(
553 Cell{pair.type + 1, pair.death}, *triangulation);
554 }
555 }
556
557 CTDiagram.reserve(pairs.size() + 1);
558
559 // find the global maximum
560 const auto nVerts = triangulation->getNumberOfVertices();
561 const SimplexId globmax = std::distance(
562 inputOffsets, std::max_element(inputOffsets, inputOffsets + nVerts));
563
564 // convert pairs to the relevant format
565 for(const auto &p : pairs) {
566 const auto isFinite = (p.death >= 0);
567 const auto death = isFinite ? p.death : globmax;
568 if(p.type == 0) {
569 const auto deathType = (isFinite && dim > 1)
572 CTDiagram.emplace_back(PersistencePair{
573 CriticalVertex{p.birth, {}, {}, {}, CriticalType::Local_minimum},
574 CriticalVertex{death, {}, {}, {}, deathType}, p.type, isFinite});
575 } else if(p.type == 1) {
576 const auto birthType
578 const auto deathType = (isFinite && dim == 3)
581 CTDiagram.emplace_back(PersistencePair{
582 CriticalVertex{p.birth, {}, {}, {}, birthType},
583 CriticalVertex{death, {}, {}, {}, deathType}, p.type, isFinite});
584 } else if(p.type == 2) {
585 CTDiagram.emplace_back(PersistencePair{
586 CriticalVertex{p.birth, {}, {}, {}, CriticalType::Saddle2},
587 CriticalVertex{death, {}, {}, {}, CriticalType::Local_maximum}, p.type,
588 isFinite});
589 }
590 }
591
592 return 0;
593}
594
595template <typename scalarType, class triangulationType>
597 std::vector<PersistencePair> &CTDiagram,
598 const scalarType *inputScalars,
599 const size_t scalarsMTime,
600 const SimplexId *inputOffsets,
601 const triangulationType *triangulation,
602 const std::vector<bool> *updateMask) {
603
604 Timer const tm{};
605 const auto dim = triangulation->getDimensionality();
606 dms_.setThreadNumber(this->threadNumber_);
607
608 dms_.buildGradient(
609 inputScalars, scalarsMTime, inputOffsets, *triangulation, updateMask);
610 std::vector<DiscreteMorseSandwich::PersistencePair> dms_pairs{};
611 dms_.computePersistencePairs(
612 dms_pairs, inputOffsets, *triangulation, this->IgnoreBoundary);
613 CTDiagram.resize(dms_pairs.size());
614
615 // transform DiscreteMorseSandwich pairs (critical cells id) to PL
616 // pairs (vertices id)
617#ifdef TTK_ENABLE_OPENMP
618#pragma omp parallel for num_threads(threadNumber_)
619#endif // TTK_ENABLE_OPENMP
620 for(size_t i = 0; i < dms_pairs.size(); ++i) {
621 auto &pair{dms_pairs[i]};
622 if(pair.type > 0) {
623 pair.birth = dms_.getCellGreaterVertex(
624 Cell{pair.type, pair.birth}, *triangulation);
625 }
626 if(pair.death != -1) {
627 pair.death = dms_.getCellGreaterVertex(
628 Cell{pair.type + 1, pair.death}, *triangulation);
629 }
630 }
631
632 // find the global maximum
633 const auto nVerts = triangulation->getNumberOfVertices();
634 const SimplexId globmax = std::distance(
635 inputOffsets, std::max_element(inputOffsets, inputOffsets + nVerts));
636
637 // convert pairs to the relevant format
638#ifdef TTK_ENABLE_OPENMP
639#pragma omp parallel for num_threads(threadNumber_)
640#endif // TTK_ENABLE_OPENMP
641 for(size_t i = 0; i < dms_pairs.size(); ++i) {
642 const auto &p{dms_pairs[i]};
643 const auto isFinite = (p.death >= 0);
644 const auto death = isFinite ? p.death : globmax;
645
646 if(p.type == 0) {
647 const auto dtype = (isFinite && dim > 1) ? CriticalType::Saddle1
649 CTDiagram[i] = PersistencePair{
650 CriticalVertex{p.birth, {}, {}, {}, CriticalType::Local_minimum},
651 CriticalVertex{death, {}, {}, {}, dtype}, p.type, isFinite};
652 } else if(p.type == 1) {
653 const auto btype
655 const auto dtype = (isFinite && dim == 3) ? CriticalType::Saddle2
657 CTDiagram[i] = PersistencePair{CriticalVertex{p.birth, {}, {}, {}, btype},
658 CriticalVertex{death, {}, {}, {}, dtype},
659 p.type, isFinite};
660 } else if(p.type == 2) {
661 const auto btype = (isFinite || dim == 3) ? CriticalType::Saddle2
663 CTDiagram[i] = PersistencePair{
664 CriticalVertex{p.birth, {}, {}, {}, btype},
665 CriticalVertex{death, {}, {}, {}, CriticalType::Local_maximum}, p.type,
666 isFinite};
667 }
668 }
669
670 return 0;
671}
672
673#if defined(TTK_ENABLE_MPI) && defined(TTK_ENABLE_OPENMP)
674template <typename scalarType, class triangulationType>
675int ttk::PersistenceDiagram::executeDiscreteMorseSandwichMPI(
676 std::vector<PersistencePair> &CTDiagram,
677 const scalarType *inputScalars,
678 const size_t scalarsMTime,
679 const SimplexId *inputOffsets,
680 const triangulationType *triangulation) {
681
682 Timer const tm{};
683 const auto dim = triangulation->getDimensionality();
684 dmsMPI_.setUseTasks(UseTasks);
685
686 dmsMPI_.buildGradient(
687 inputScalars, scalarsMTime, inputOffsets, *triangulation);
688 std::vector<DiscreteMorseSandwichMPI::PersistencePair> dms_pairs{};
689 dmsMPI_.computePersistencePairs(
690 dms_pairs, inputOffsets, *triangulation, this->IgnoreBoundary);
691 CTDiagram.resize(dms_pairs.size());
692
693 // transform DiscreteMorseSandwich pairs (critical cells id) to PL
694 // pairs (vertices id)
695 struct dataRequest {
696 ttk::SimplexId gid_;
697 ttk::SimplexId lid_;
698 char dim_;
699 char isBirth_;
700 };
701 // find the global maximum
702 const auto nVerts = triangulation->getNumberOfVertices();
703 const SimplexId localMaxId = std::distance(
704 inputOffsets, std::max_element(inputOffsets, inputOffsets + nVerts));
705 MPI_Datatype MPI_SimplexId = getMPIType(localMaxId);
706 struct {
707 long int offset;
708 int rank;
709 } localMax, globalMax;
710 localMax.rank = triangulation->getVertexRank(localMaxId);
711 localMax.offset = static_cast<long int>(inputOffsets[localMaxId]);
712 MPI_Allreduce(
713 &localMax, &globalMax, 1, MPI_LONG_INT, MPI_MAXLOC, ttk::MPIcomm_);
714 ttk::SimplexId globmax{-1};
715 double maxMetaData[4];
716 if(globalMax.rank == ttk::MPIrank_) {
717 float coords[3];
718 globmax = triangulation->getVertexGlobalId(localMaxId);
719 triangulation->getVertexPoint(localMaxId, coords[0], coords[1], coords[2]);
720 maxMetaData[0] = coords[0];
721 maxMetaData[1] = coords[1];
722 maxMetaData[2] = coords[2];
723 maxMetaData[3] = inputScalars[localMaxId];
724 }
725 MPI_Bcast(&globmax, 1, MPI_SimplexId, globalMax.rank, ttk::MPIcomm_);
726 MPI_Bcast(maxMetaData, 4, MPI_DOUBLE, globalMax.rank, ttk::MPIcomm_);
727
728 std::vector<std::vector<dataRequest>> sendRecvBuffer(ttk::MPIsize_);
729
730 const auto createDataRequestMPIType =
731 [&MPI_SimplexId](MPI_Datatype &MPI_MessageType) {
732 MPI_Datatype types[] = {MPI_SimplexId, MPI_SimplexId, MPI_CHAR, MPI_CHAR};
733 int lengths[] = {1, 1, 1, 1};
734 const long int mpi_offsets[]
735 = {offsetof(dataRequest, gid_), offsetof(dataRequest, lid_),
736 offsetof(dataRequest, dim_), offsetof(dataRequest, isBirth_)};
737 MPI_Type_create_struct(4, lengths, mpi_offsets, types, &MPI_MessageType);
738 MPI_Type_commit(&MPI_MessageType);
739 };
740
741 struct dataResponse {
742 ttk::SimplexId lid_{-1};
743 ttk::SimplexId vertexGid_{-1};
744 ttk::SimplexId offset_{-1};
745 float coords_[3] = {0, 0, 0};
746 double sfValue_{0};
747 char isBirth_{0};
748 };
749
750 const auto createDataResponseMPIType =
751 [&MPI_SimplexId](MPI_Datatype &MPI_MessageType) {
752 MPI_Datatype types[] = {MPI_SimplexId, MPI_SimplexId, MPI_SimplexId,
753 MPI_FLOAT, MPI_DOUBLE, MPI_CHAR};
754 int lengths[] = {1, 1, 1, 3, 1, 1};
755 const long int mpi_offsets[]
756 = {offsetof(dataResponse, lid_), offsetof(dataResponse, vertexGid_),
757 offsetof(dataResponse, offset_), offsetof(dataResponse, coords_),
758 offsetof(dataResponse, sfValue_), offsetof(dataResponse, isBirth_)};
759 MPI_Type_create_struct(6, lengths, mpi_offsets, types, &MPI_MessageType);
760 MPI_Type_commit(&MPI_MessageType);
761 };
762
763 const auto fillBirthData
764 = [&dim](PersistencePair &CTPair,
765 DiscreteMorseSandwichMPI::PersistencePair &p,
766 ttk::SimplexId birthId) {
767 CTPair.birth.id = birthId;
768 if(p.type == 0) {
769 CTPair.birth.type = CriticalType::Local_minimum;
770 } else if(p.type == 1) {
771 CTPair.birth.type
773 } else if(p.type == 2) {
774 CTPair.birth.type = ((p.death >= 0) || dim == 3)
777 }
778 };
779
780 const auto augmentBirthPersistence =
781 [&triangulation, &inputOffsets](
782 PersistencePair &CTPair, ttk::SimplexId lid, const scalarType *scalars) {
783 triangulation->getVertexPoint(lid, CTPair.birth.coords[0],
784 CTPair.birth.coords[1],
785 CTPair.birth.coords[2]);
786 CTPair.birth.sfValue = scalars[lid];
787 CTPair.birth.offset = inputOffsets[lid];
788 };
789
790 const auto fillDeathData = [&dim](
791 PersistencePair &CTPair,
792 DiscreteMorseSandwichMPI::PersistencePair &p,
793 ttk::SimplexId deathId) {
794 const auto isFinite = (p.death >= 0);
795 CTPair.death.id = deathId;
796 if(p.type == 0) {
797 CTPair.death.type = (isFinite && dim > 1) ? CriticalType::Saddle1
799 } else if(p.type == 1) {
800 CTPair.death.type = (isFinite && dim == 3) ? CriticalType::Saddle2
802 } else if(p.type == 2) {
803 CTPair.death.type = CriticalType::Local_maximum;
804 }
805 };
806
807 const auto augmentDeathPersistence =
808 [&triangulation, &inputOffsets](
809 PersistencePair &CTPair, ttk::SimplexId lid, const scalarType *scalars) {
810 triangulation->getVertexPoint(lid, CTPair.death.coords[0],
811 CTPair.death.coords[1],
812 CTPair.death.coords[2]);
813 CTPair.death.sfValue = scalars[lid];
814 CTPair.death.offset = inputOffsets[lid];
815 };
816
817 const auto getBirthSimplexType = [&dim](const int type) {
818 switch(type) {
819 case 0:
820 return 0;
821 case 1:
822 return 1;
823 default:
824 if(dim == 3) {
825 return 2;
826 } else {
827 return 1;
828 }
829 }
830 };
831
832 const auto getDeathSimplexType = [&dim](const int type) {
833 switch(type) {
834 case 0:
835 return 1;
836 case 1:
837 return 2;
838 default:
839 if(dim == 3) {
840 return 3;
841 } else {
842 return 2;
843 }
844 }
845 };
846 for(ttk::SimplexId i = 0; i < static_cast<ttk::SimplexId>(dms_pairs.size());
847 ++i) {
848 auto &pair{dms_pairs[i]};
849 int simplexType = getBirthSimplexType(pair.type);
851 = triangulation->getSimplexLocalId(pair.birth, simplexType);
852 if(lid != -1
853 && triangulation->getSimplexRank(lid, simplexType) == ttk::MPIrank_) {
854 if(pair.type > 0) {
855 lid
856 = dmsMPI_.getCellGreaterVertex(Cell{pair.type, lid}, *triangulation);
857 pair.birth = triangulation->getVertexGlobalId(lid);
858 }
859 CTDiagram[i].dim = pair.type;
860 CTDiagram[i].isFinite = (pair.death >= 0);
861 // Add all the other stuff
862 fillBirthData(CTDiagram[i], pair, pair.birth);
863 augmentBirthPersistence(CTDiagram[i], lid, inputScalars);
864 } else {
865 sendRecvBuffer[ttk::MPIrank_].emplace_back(
866 dataRequest{pair.birth, i, static_cast<char>(pair.type), 1});
867 }
868 if(pair.death == -1) {
869 CTDiagram[i].dim = pair.type;
870 CTDiagram[i].isFinite = (pair.death >= 0);
871 pair.death = globmax;
872 fillDeathData(CTDiagram[i], pair, pair.death);
873 CTDiagram[i].death.coords[0] = maxMetaData[0];
874 CTDiagram[i].death.coords[1] = maxMetaData[1];
875 CTDiagram[i].death.coords[2] = maxMetaData[2];
876 CTDiagram[i].death.sfValue = maxMetaData[3];
877 CTDiagram[i].death.offset = globalMax.offset;
878 } else {
879 simplexType = getDeathSimplexType(pair.type);
880 lid = triangulation->getSimplexLocalId(pair.death, simplexType);
881 if(lid != -1
882 && triangulation->getSimplexRank(lid, simplexType) == ttk::MPIrank_) {
883 CTDiagram[i].dim = pair.type;
884 CTDiagram[i].isFinite = (pair.death >= 0);
885 lid = dmsMPI_.getCellGreaterVertex(
886 Cell{pair.type + 1, lid}, *triangulation);
887 pair.death = triangulation->getVertexGlobalId(lid);
888 fillDeathData(CTDiagram[i], pair, pair.death);
889 augmentDeathPersistence(CTDiagram[i], lid, inputScalars);
890 } else {
891 sendRecvBuffer[ttk::MPIrank_].emplace_back(
892 dataRequest{pair.death, i, static_cast<char>(pair.type), 0});
893 }
894 }
895 }
896 // Broadcast all the data to everyone
897 // First, broadcast the size of the data to send
898 std::vector<int> recvBufferSize(ttk::MPIsize_, 0);
899 recvBufferSize[ttk::MPIrank_] = sendRecvBuffer[ttk::MPIrank_].size();
900 for(int i = 0; i < ttk::MPIsize_; i++) {
901 MPI_Bcast(&recvBufferSize[i], 1, MPI_INTEGER, i, ttk::MPIcomm_);
902 if(i != ttk::MPIrank_) {
903 sendRecvBuffer[i].resize(recvBufferSize[i]);
904 }
905 }
906 // Then, broadcast the actual data
907 MPI_Datatype MPI_requestDataType;
908 createDataRequestMPIType(MPI_requestDataType);
909 for(int i = 0; i < ttk::MPIsize_; i++) {
910 MPI_Bcast(sendRecvBuffer[i].data(), recvBufferSize[i], MPI_requestDataType,
911 i, ttk::MPIcomm_);
912 }
913 std::vector<std::vector<dataResponse>> response(ttk::MPIsize_);
914
915 // For each received element:
916 // if it is own by the triangulation, get the data and place to send back
917 for(int i = 0; i < ttk::MPIsize_; i++) {
918 if(i != ttk::MPIrank_) {
919 for(int j = 0; j < recvBufferSize[i]; j++) {
920 auto &element{sendRecvBuffer[i][j]};
921 int simplexType;
922 if(element.isBirth_) {
923 simplexType = getBirthSimplexType(element.dim_);
924 } else {
925 simplexType = getDeathSimplexType(element.dim_);
926 }
928 = triangulation->getSimplexLocalId(element.gid_, simplexType);
929
930 if(lid != -1
931 && triangulation->getSimplexRank(lid, simplexType)
932 == ttk::MPIrank_) {
933 // Add the relevant data
934 struct dataResponse res {
935 .lid_ = element.lid_, .isBirth_ = element.isBirth_
936 };
937 ttk::SimplexId vLid = dmsMPI_.getCellGreaterVertex(
938 Cell{element.dim_ + (1 - element.isBirth_), lid}, *triangulation);
939 res.vertexGid_ = triangulation->getVertexGlobalId(vLid);
940 res.offset_ = inputOffsets[vLid];
941 triangulation->getVertexPoint(
942 vLid, res.coords_[0], res.coords_[1], res.coords_[2]);
943 res.sfValue_ = inputScalars[vLid];
944 // Store to send
945 response[i].emplace_back(res);
946 }
947 }
948 }
949 }
950 // Send back the data
951 // First, exchange the size of the data to exchange
952 std::vector<dataResponse> responseBuffer;
953 std::vector<int> sendBufferSize(ttk::MPIsize_, 0);
954 std::vector<int> sendDispls(ttk::MPIsize_, 0);
955 std::vector<int> recvDispls(ttk::MPIsize_, 0);
956 std::vector<dataResponse> recvBuffer;
957
958 sendBufferSize[0] = response[0].size();
959 responseBuffer.insert(
960 responseBuffer.end(), response[0].begin(), response[0].end());
961 for(int i = 1; i < ttk::MPIsize_; i++) {
962 sendBufferSize[i] = response[i].size();
963 sendDispls[i] = sendDispls[i - 1] + sendBufferSize[i - 1];
964 responseBuffer.insert(
965 responseBuffer.end(), response[i].begin(), response[i].end());
966 }
967 MPI_Alltoall(sendBufferSize.data(), 1, MPI_INTEGER, recvBufferSize.data(), 1,
968 MPI_INTEGER, ttk::MPIcomm_);
969 for(int i = 1; i < ttk::MPIsize_; i++) {
970 recvDispls[i] = recvDispls[i - 1] + recvBufferSize[i - 1];
971 }
972 recvBuffer.resize(recvDispls.back() + recvBufferSize.back());
973 MPI_Datatype MPI_responseDataType;
974 createDataResponseMPIType(MPI_responseDataType);
975
976 // Then, exchange to actual data
977 MPI_Alltoallv(responseBuffer.data(), sendBufferSize.data(), sendDispls.data(),
978 MPI_responseDataType, recvBuffer.data(), recvBufferSize.data(),
979 recvDispls.data(), MPI_responseDataType, ttk::MPIcomm_);
980 // Receive the data and store it appropriately
981 for(const auto &element : recvBuffer) {
982 auto &CTPair{CTDiagram[element.lid_]};
983 if(element.isBirth_) {
984 fillBirthData(CTPair, dms_pairs[element.lid_], element.vertexGid_);
985 CTPair.birth.coords[0] = element.coords_[0];
986 CTPair.birth.coords[1] = element.coords_[1];
987 CTPair.birth.coords[2] = element.coords_[2];
988 CTPair.birth.sfValue = element.sfValue_;
989 CTPair.birth.offset = element.offset_;
990 } else {
991 fillDeathData(CTPair, dms_pairs[element.lid_], element.vertexGid_);
992 CTPair.death.coords[0] = element.coords_[0];
993 CTPair.death.coords[1] = element.coords_[1];
994 CTPair.death.coords[2] = element.coords_[2];
995 CTPair.death.sfValue = element.sfValue_;
996 CTPair.death.offset = element.offset_;
997 }
998 }
999 // Create MPI type for Critical Vertex
1000 const auto createCriticalVertexMPIType =
1001 [&MPI_SimplexId](MPI_Datatype &MPI_MessageType) {
1002 MPI_Datatype types[]
1003 = {MPI_SimplexId, MPI_DOUBLE, MPI_SimplexId, MPI_FLOAT, MPI_INTEGER};
1004 int lengths[] = {1, 1, 1, 3, 1};
1005 const long int mpi_offsets[]
1006 = {offsetof(CriticalVertex, id), offsetof(CriticalVertex, sfValue),
1007 offsetof(CriticalVertex, offset), offsetof(CriticalVertex, coords),
1008 offsetof(CriticalVertex, type)};
1009 MPI_Type_create_struct(5, lengths, mpi_offsets, types, &MPI_MessageType);
1010 // printMsg("create_struct done");
1011 MPI_Type_commit(&MPI_MessageType);
1012 };
1013 MPI_Datatype MPI_CriticalVertex;
1014 createCriticalVertexMPIType(MPI_CriticalVertex);
1015 // Create MPI type for PersistencePair
1016 const auto createPersistencePairMPIType =
1017 [&MPI_CriticalVertex, &MPI_SimplexId](MPI_Datatype &MPI_MessageType) {
1018 MPI_Datatype types[]
1019 = {MPI_CriticalVertex, MPI_CriticalVertex, MPI_SimplexId, MPI_CHAR};
1020 int lengths[] = {1, 1, 1, 1};
1021 const long int mpi_offsets[]
1022 = {offsetof(PersistencePair, birth), offsetof(PersistencePair, death),
1023 offsetof(PersistencePair, dim), offsetof(PersistencePair, isFinite)};
1024 MPI_Type_create_struct(4, lengths, mpi_offsets, types, &MPI_MessageType);
1025 MPI_Type_commit(&MPI_MessageType);
1026 };
1027 MPI_Datatype MPI_PersistencePair;
1028 createPersistencePairMPIType(MPI_PersistencePair);
1029 // Sort in parallel using the offset of the birth and psort
1030 std::vector<ttk::SimplexId> vertexDistribution(ttk::MPIsize_);
1031 ttk::SimplexId localVertexNumber = CTDiagram.size();
1032 MPI_Allgather(&localVertexNumber, 1, MPI_SimplexId, vertexDistribution.data(),
1033 1, MPI_SimplexId, ttk::MPIcomm_);
1034 p_sort::parallel_sort<PersistencePair>(
1036 vertexDistribution, MPI_PersistencePair, MPI_SimplexId, threadNumber_);
1037 return 0;
1038};
1039#endif
1040
1041template <typename scalarType, class triangulationType>
1043 std::vector<PersistencePair> &CTDiagram,
1044 const scalarType *inputScalars,
1045 const triangulationType *triangulation) {
1046
1047 approxT_.setDebugLevel(debugLevel_);
1048 approxT_.setThreadNumber(threadNumber_);
1049 approxT_.setupTriangulation(const_cast<ttk::ImplicitTriangulation *>(
1050 (const ImplicitTriangulation *)triangulation));
1051 approxT_.setStartingResolutionLevel(StartingResolutionLevel);
1052 approxT_.setStoppingResolutionLevel(StoppingResolutionLevel);
1053 approxT_.setPreallocateMemory(true);
1054 approxT_.setEpsilon(Epsilon);
1055
1056 std::vector<ApproximateTopology::PersistencePair> resultDiagram{};
1057
1058 approxT_.computeApproximatePD(
1059 resultDiagram, inputScalars, (scalarType *)outputScalars_,
1061
1062 // create the final diagram
1063 for(const auto &p : resultDiagram) {
1064 if(p.pairType == 0) {
1065 CTDiagram.emplace_back(PersistencePair{
1066 CriticalVertex{p.birth, {}, {}, {}, CriticalType::Local_minimum},
1067 CriticalVertex{p.death, {}, {}, {}, CriticalType::Saddle1}, p.pairType,
1068 true});
1069 } else if(p.pairType == 2) {
1070 CTDiagram.emplace_back(PersistencePair{
1071 CriticalVertex{p.birth, {}, {}, {}, CriticalType::Saddle2},
1072 CriticalVertex{p.death, {}, {}, {}, CriticalType::Local_maximum},
1073 p.pairType, true});
1074 } else if(p.pairType == -1) {
1075 CTDiagram.emplace_back(PersistencePair{
1076 CriticalVertex{p.birth, {}, {}, {}, CriticalType::Local_minimum},
1077 CriticalVertex{p.death, {}, {}, {}, CriticalType::Local_maximum},
1078 p.pairType, false});
1079 }
1080 }
1081
1082 return 0;
1083}
1084
1085template <class triangulationType>
1087 std::vector<PersistencePair> &CTDiagram,
1088 const SimplexId *inputOffsets,
1089 const triangulationType *triangulation) {
1090
1091 progT_.setDebugLevel(debugLevel_);
1092 progT_.setThreadNumber(threadNumber_);
1093 progT_.setupTriangulation(const_cast<ttk::ImplicitTriangulation *>(
1094 (const ImplicitTriangulation *)triangulation));
1095 progT_.setStartingResolutionLevel(StartingResolutionLevel);
1096 progT_.setStoppingResolutionLevel(StoppingResolutionLevel);
1097 progT_.setTimeLimit(TimeLimit);
1098 progT_.setIsResumable(IsResumable);
1099 progT_.setPreallocateMemory(true);
1100
1101 std::vector<ProgressiveTopology::PersistencePair> resultDiagram{};
1102
1103 progT_.computeProgressivePD(resultDiagram, inputOffsets);
1104
1105 // create the final diagram
1106 for(const auto &p : resultDiagram) {
1107 if(p.pairType == 0) {
1108 CTDiagram.emplace_back(PersistencePair{
1109 CriticalVertex{p.birth, {}, {}, {}, CriticalType::Local_minimum},
1110 CriticalVertex{p.death, {}, {}, {}, CriticalType::Saddle1}, p.pairType,
1111 true});
1112 } else if(p.pairType == 2) {
1113 CTDiagram.emplace_back(PersistencePair{
1114 CriticalVertex{p.birth, {}, {}, {}, CriticalType::Saddle2},
1115 CriticalVertex{p.death, {}, {}, {}, CriticalType::Local_maximum},
1116 p.pairType, true});
1117 } else if(p.pairType == -1) {
1118 CTDiagram.emplace_back(PersistencePair{
1119 CriticalVertex{p.birth, {}, {}, {}, CriticalType::Local_minimum},
1120 CriticalVertex{p.death, {}, {}, {}, CriticalType::Local_maximum}, 0,
1121 false});
1122 }
1123 }
1124
1125 return 0;
1126}
1127
1128template <typename scalarType, class triangulationType>
1130 std::vector<PersistencePair> &CTDiagram,
1131 const scalarType *inputScalars,
1132 const SimplexId *inputOffsets,
1133 const triangulationType *triangulation) {
1134
1135 contourTree_.setVertexScalars(inputScalars);
1137 contourTree_.setVertexSoSoffsets(inputOffsets);
1138 contourTree_.setSegmentation(false);
1139 contourTree_.build<scalarType>(triangulation);
1140
1141 // get persistence pairs
1142 std::vector<std::tuple<ttk::SimplexId, ttk::SimplexId, scalarType>> JTPairs;
1143 std::vector<std::tuple<ttk::SimplexId, ttk::SimplexId, scalarType>> STPairs;
1144 contourTree_.computePersistencePairs<scalarType>(JTPairs, true);
1145 contourTree_.computePersistencePairs<scalarType>(STPairs, false);
1146
1147 // merge pairs
1148 const auto JTSize = JTPairs.size();
1149 const auto STSize = STPairs.size();
1150 std::vector<std::tuple<ttk::SimplexId, ttk::SimplexId, scalarType, bool>>
1151 CTPairs(JTSize + STSize);
1152 for(size_t i = 0; i < JTSize; ++i) {
1153 const auto &x = JTPairs[i];
1154 CTPairs[i]
1155 = std::make_tuple(std::get<0>(x), std::get<1>(x), std::get<2>(x), true);
1156 }
1157 for(size_t i = 0; i < STSize; ++i) {
1158 const auto &x = STPairs[i];
1159 CTPairs[JTSize + i]
1160 = std::make_tuple(std::get<0>(x), std::get<1>(x), std::get<2>(x), false);
1161 }
1162
1163 // remove the last pair which is present two times (global extrema pair)
1164 if(!CTPairs.empty()) {
1165 auto cmp =
1166 [](
1167 const std::tuple<ttk::SimplexId, ttk::SimplexId, scalarType, bool> &a,
1168 const std::tuple<ttk::SimplexId, ttk::SimplexId, scalarType, bool> &b) {
1169 return std::get<2>(a) < std::get<2>(b);
1170 };
1171
1172 std::sort(CTPairs.begin(), CTPairs.end(), cmp);
1173 CTPairs.erase(CTPairs.end() - 1);
1174 }
1175
1176 // get persistence diagrams
1178
1179 return 0;
1180}
1181
1182template <class triangulationType>
1184 const triangulationType *ttkNotUsed(triangulation)) {
1185
1188 && !std::is_same<ttk::ImplicitWithPreconditions, triangulationType>::value
1189 && !std::is_same<ttk::ImplicitNoPreconditions, triangulationType>::value) {
1190
1191 printWrn("Explicit, Compact or Periodic triangulation detected.");
1192 printWrn("Defaulting to the FTM backend.");
1193
1195 }
1196}
1197
1198template <class triangulationType>
1200 const triangulationType *const triangulation) {
1201
1204 return;
1205 }
1206
1207 if(!triangulation->isManifold()) {
1208 this->printWrn("Non-manifold data-set detected.");
1209 this->printWrn("Defaulting to the Persistence Simplex backend.");
1210
1212 }
1213}
#define ttkNotUsed(x)
Mark function/method parameters that are not used in the function body at all.
Definition BaseClass.h:47
AbstractTriangulation is an interface class that defines an interface for efficient traversal methods...
TTK processing package for progressive Topological Data Analysis.
int debugLevel_
Definition Debug.h:379
int printWrn(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition Debug.h:159
int printErr(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition Debug.h:149
TTK DiscreteMorseSandwichMPI processing package.
TTK DiscreteMorseSandwich processing package.
ImplicitTriangulation is a class that provides time and memory efficient traversal methods on triangu...
int executePersistentSimplex(std::vector< PersistencePair > &CTDiagram, const SimplexId *inputOffsets, const triangulationType *triangulation)
void setComputeSadSad(const bool data)
void setOutputMonotonyOffsets(void *data)
void sortPersistenceDiagram(std::vector< PersistencePair > &diagram, const SimplexId *const offsets) const
void preconditionTriangulation(AbstractTriangulation *triangulation)
ttk::ProgressiveTopology progT_
void setDeltaApproximate(double data)
void setBackend(const BACKEND be)
ttk::CriticalType getNodeType(ftm::FTMTree_MT *tree, ftm::TreeType treeType, const SimplexId vertexId) const
void checkProgressivityRequirement(const triangulationType *triangulation)
void setOutputScalars(void *data)
int executeFTM(std::vector< PersistencePair > &CTDiagram, const scalarType *inputScalars, const SimplexId *inputOffsets, const triangulationType *triangulation)
dcg::DiscreteGradient dcg_
int executeProgressiveTopology(std::vector< PersistencePair > &CTDiagram, const SimplexId *inputOffsets, const triangulationType *triangulation)
int execute(std::vector< PersistencePair > &CTDiagram, const scalarType *inputScalars, const size_t scalarsMTime, const SimplexId *inputOffsets, const triangulationType *triangulation, const std::vector< bool > *updateMask=nullptr)
ttk::ApproximateTopology approxT_
DiscreteMorseSandwich dms_
int executeDiscreteMorseSandwich(std::vector< PersistencePair > &CTDiagram, const scalarType *inputScalars, const size_t scalarsMTime, const SimplexId *inputOffsets, const triangulationType *triangulation, const std::vector< bool > *updateMask=nullptr)
void setComputeMinSad(const bool data)
void setOutputOffsets(void *data)
void augmentPersistenceDiagram(std::vector< PersistencePair > &persistencePairs, const scalarType *const scalars, const triangulationType *triangulation)
Complete a ttk::DiagramType instance with scalar field values (useful for persistence) and 3D coordin...
void setComputeSadMax(const bool data)
int executeApproximateTopology(std::vector< PersistencePair > &CTDiagram, const scalarType *inputScalars, const triangulationType *triangulation)
PersistentSimplexPairs psp_
void checkManifold(const triangulationType *const triangulation)
int computeCTPersistenceDiagram(ftm::FTMTreePP &tree, const std::vector< std::tuple< ttk::SimplexId, ttk::SimplexId, scalarType, bool > > &pairs, std::vector< PersistencePair > &diagram) const
Textbook algorithm to find persistence pairs.
TTK processing package for progressive Topological Data Analysis.
TTK discreteGradient processing package.
FTMTree_MT * getSplitTree()
Definition FTMTree_CT.h:51
FTMTree_MT * getJoinTree()
Definition FTMTree_CT.h:47
bool comp(const PersistencePair a, const PersistencePair b)
bool oppositeComp(const PersistencePair a, const PersistencePair b)
TTK base package defining the standard types.
COMMON_EXPORTS int MPIsize_
Definition BaseClass.cpp:10
int SimplexId
Identifier type for simplices of any dimension.
Definition DataTypes.h:22
CriticalType
default value for critical index
Definition DataTypes.h:88
COMMON_EXPORTS int MPIrank_
Definition BaseClass.cpp:9
printMsg(debug::output::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/| (_) |"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)