Merge pull request #1112 from eve-le-guillou/mpi_DMS_pr

[MPI] Distributed Discrete Morse Sandwich: Hybrid MPI+thread Persistence Diagram Computation

Author: julien-tierny

.github/workflows/package.yml

@@ -164,7 +164,7 @@ jobs:
         rm python/nestedTrackingFromOverlap.py
         # remove examples which fill up the memory
         rm python/topologicalOptimization_darkSky.py
-        python3 -u python/run.py
+        pvpython -u python/run.py
 
     - name: Test ttk-data Python scripts results [NOT ENFORCED]
       continue-on-error: true

.github/workflows/test.yml

@@ -247,7 +247,7 @@ jobs:
         rm python/mpiExample.py
         # remove examples which fill up the memory
         rm python/topologicalOptimization_darkSky.py
-        python3 -u python/run.py
+        pvpython -u python/run.py
 
     - name: Test ttk-data Python scripts results [NOT ENFORCED]
       if: matrix.testSet == 'pyTests'
@@ -398,6 +398,7 @@ jobs:
         rm python/nestedTrackingFromOverlap.py
         # tmp removal of mpi example
         rm python/mpiExample.py
+        rm python/distributedPersistenceDiagram.py
         # remove examples which fill up the memory
         rm python/topologicalOptimization_darkSky.py
         # remove sklearn examples -- buggy on MacOS since sklearn-1.7

CMake/config.cmake

@@ -46,7 +46,7 @@ if (TTK_ENABLE_MPI)
   option(TTK_ENABLE_MPI_TIME "Enable time measuring for MPI computation" FALSE)
   mark_as_advanced(TTK_ENABLE_MPI_TIME)
   option(TTK_ENABLE_MPI_RANK_ID_INT "Enable rank ids of type int (default char) for distributed sort" FALSE)
-  mark_as_advanced(TTK_ENABLE_MPI_RANK_ID_TIME)
+  mark_as_advanced(TTK_ENABLE_MPI_RANK_ID_INT)
 
 endif()
 

ChangeLog.md

@@ -1,7 +1,9 @@
 ## TTK - ChangeLog
 =
 ### dev
+- Distributed computation of persistent homology!
 - New backend for TrackingFromFields (critical point based)
+- Fast planar Rips filtration persistence computation
 
 ### 1.3.0
 - Wasserstein Auto-Encoders of Merge Trees / Pers. Diagrams (IEEE TVCG 2024)

core/base/abstractTriangulation/AbstractTriangulation.h

@@ -2828,6 +2828,23 @@ namespace ttk {
       return -1;
     }
 
+    inline ttk::SimplexId getSimplexLocalId(const SimplexId gsid,
+                                            const int type) const {
+      if(!ttk::isRunningWithMPI()) {
+        return gsid;
+      }
+      switch(type) {
+        case 0:
+          return this->getVertexLocalId(gsid);
+        case 1:
+          return this->getEdgeLocalId(gsid);
+        case 2:
+          return this->getTriangleLocalId(gsid);
+        default:
+          return this->getCellLocalId(gsid);
+      }
+    }
+
     virtual inline int getVertexRank(const SimplexId lvid) const {
 
       if(!ttk::isRunningWithMPI()) {
@@ -2852,6 +2869,54 @@ namespace ttk {
       return this->getVertexRankInternal(lvid);
     }
 
+    virtual inline int getEdgeRank(const SimplexId lvid) const {
+
+      if(!ttk::isRunningWithMPI()) {
+        return 0;
+      }
+#ifndef TTK_ENABLE_KAMIKAZE
+      if(this->getDimensionality() != 1 && this->getDimensionality() != 2
+         && this->getDimensionality() != 3) {
+        this->printErr("Only 1D, 2D and 3D datasets are supported");
+        return -1;
+      }
+      if(!this->hasPreconditionedDistributedEdges_) {
+        this->printErr("EdgeRankId query without pre-process!");
+        this->printErr(
+          "Please call preconditionDistributedEdges() in a pre-process.");
+        return -1;
+      }
+      if(lvid < 0 || lvid >= this->getNumberOfEdges()) {
+        return -1;
+      }
+#endif // TTK_ENABLE_KAMIKAZE
+      return this->getEdgeRankInternal(lvid);
+    }
+
+    virtual inline int getTriangleRank(const SimplexId lvid) const {
+
+      if(!ttk::isRunningWithMPI()) {
+        return 0;
+      }
+#ifndef TTK_ENABLE_KAMIKAZE
+      if(this->getDimensionality() != 1 && this->getDimensionality() != 2
+         && this->getDimensionality() != 3) {
+        this->printErr("Only 1D, 2D and 3D datasets are supported");
+        return -1;
+      }
+      if(!this->hasPreconditionedDistributedTriangles_) {
+        this->printErr("TriangleRankId query without pre-process!");
+        this->printErr(
+          "Please call preconditionDistributedTriangles() in a pre-process.");
+        return -1;
+      }
+      if(lvid < 0 || lvid >= this->getNumberOfTriangles()) {
+        return -1;
+      }
+#endif // TTK_ENABLE_KAMIKAZE
+      return this->getTriangleRankInternal(lvid);
+    }
+
     virtual inline int getCellRank(const SimplexId lcid) const {
 #ifndef TTK_ENABLE_KAMIKAZE
       if(this->getDimensionality() != 1 && this->getDimensionality() != 2
@@ -2870,11 +2935,31 @@ namespace ttk {
       }
 #endif // TTK_ENABLE_KAMIKAZE
       if(!ttk::isRunningWithMPI()) {
-        return lcid;
+        return ttk::MPIrank_;
       }
       return this->getCellRankInternal(lcid);
     }
 
+    inline int getSimplexRank(const SimplexId lsid, const int type) const {
+      if(!ttk::isRunningWithMPI()) {
+        return ttk::MPIrank_;
+      }
+      switch(type) {
+        case 0:
+          return this->getVertexRank(lsid);
+        case 1:
+          return this->getEdgeRank(lsid);
+        case 2: {
+          if(getDimensionality() == 3)
+            return this->getTriangleRank(lsid);
+          return this->getCellRank(lsid);
+        }
+        case 3:
+          return this->getCellRank(lsid);
+      }
+      return -1;
+    }
+
     virtual inline const std::vector<int> &getNeighborRanks() const {
       return this->neighborRanks_;
     }
@@ -3046,6 +3131,16 @@ namespace ttk {
       return -1;
     }
 
+    virtual inline int
+      getEdgeRankInternal(const SimplexId ttkNotUsed(lvid)) const {
+      return -1;
+    }
+
+    virtual inline int
+      getTriangleRankInternal(const SimplexId ttkNotUsed(lvid)) const {
+      return -1;
+    }
+
     virtual inline int
       getCellRankInternal(const SimplexId ttkNotUsed(lcid)) const {
       return -1;

core/base/arrayPreconditioning/ArrayPreconditioning.h

@@ -344,9 +344,10 @@ namespace ttk {
 
         ttk::SimplexId verticesToSortSize = verticesToSort.size();
         // Compute the order of the first element of the current process
-        ttk::SimplexId orderOffset
-          = std::accumulate(vertexDistribution.begin(),
-                            vertexDistribution.begin() + ttk::MPIrank_, 0);
+        ttk::SimplexId orderOffset{0};
+        for(int i = 0; i < ttk::MPIrank_; i++) {
+          orderOffset += vertexDistribution[i];
+        }
         // nbChunk, rest and nbChunkTotal are used to compute the
         // post-processing bit by bit.
         ttk::SimplexId nbChunk = std::floor(verticesToSortSize / ChunkSize);
@@ -412,7 +413,7 @@ namespace ttk {
     }
 
   protected:
-    bool GlobalOrder{false};
+    bool GlobalOrder{true};
     // This value has been chosen for systems of 128 Gb of memory per computing
     // node. For systems with much smaller memory, it may be inadequate and
     // require a smaller value.