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