Merge pull request #1105 from MClemot/TopologicallyConstrainedDR

[DimensionReduction] New backend

Author: julien-tierny

ChangeLog.md

@@ -1,7 +1,8 @@
 ## TTK - ChangeLog
 =
 ### dev
-- Distributed computation of persistent homology!
+- Cycle-aware dimensionality reduction (TopoAE++)
+- Distributed computation of persistent homology! (IEEE TPDS 2025)
 - New backend for TrackingFromFields (critical point based)
 - Fast planar Rips filtration persistence computation
 - Migration to ParaView 6

core/base/dimensionReduction/CMakeLists.txt

@@ -6,6 +6,7 @@ ttk_add_base_library(dimensionReduction
   DEPENDS
     triangulation
     topoMap
+    topologicalDimensionReduction
   )
 
 install(

core/base/dimensionReduction/DimensionReduction.cpp

@@ -1,8 +1,6 @@
 #include <DimensionReduction.h>
 #include <TopoMap.h>
 
-#include <map>
-
 #define VALUE_TO_STRING(x) #x
 #define VALUE(x) VALUE_TO_STRING(x)
 
@@ -77,6 +75,40 @@ int DimensionReduction::execute(
     return 0;
   }
 
+  if(this->Method == METHOD::AE) {
+#ifdef TTK_ENABLE_TORCH
+    TopologicalDimensionReduction tcdr(
+      ae_CUDA, ae_Deterministic, ae_Seed, NumberOfComponents, ae_Epochs,
+      ae_LearningRate, ae_Optimizer, ae_Method, ae_Model, ae_Architecture,
+      ae_Activation, ae_BatchSize, ae_BatchNormalization, ae_RegCoefficient,
+      IsInputImages);
+    tcdr.setDebugLevel(debugLevel_);
+    tcdr.setThreadNumber(threadNumber_);
+
+    outputEmbedding.resize(NumberOfComponents);
+    for(int d = 0; d < NumberOfComponents; d++)
+      outputEmbedding[d].resize(nRows);
+
+    if(ae_PreOptimize) {
+      DimensionReduction initDR;
+      initDR.setDebugLevel(debugLevel_);
+      initDR.setInputMethod(ae_PreOptimizeMethod);
+      std::vector<std::vector<double>> latentInitialization;
+      initDR.execute(latentInitialization, inputMatrix, nRows, nColumns);
+      tcdr.setLatentInitialization(latentInitialization);
+    }
+
+    tcdr.execute(outputEmbedding, inputMatrix, nRows);
+
+    this->printMsg("Computed AE dimension reduction", 1.0, t.getElapsedTime(),
+                   threadNumber_);
+    return 0;
+#else
+    this->printErr("Unavailable backend: Torch is required.");
+    return 1;
+#endif
+  }
+
 #ifdef TTK_ENABLE_SCIKIT_LEARN
 #ifndef TTK_ENABLE_KAMIKAZE
   if(majorVersion_ < '3')

core/base/dimensionReduction/DimensionReduction.h

@@ -40,6 +40,9 @@
 ///   href="https://topology-tool-kit.github.io/examples/persistentGenerators_periodicPicture/">Persistent
 ///   Generators Periodic Picture example</a> \n
 ///   - <a
+///   href="https://topology-tool-kit.github.io/examples/topoAEppTeaser/">Topological
+///   Autoencoders++ Teaser example</a> \n
+///   - <a
 ///   href="https://topology-tool-kit.github.io/examples/topoMapTeaser/">TopoMap
 ///   Teaser example</a> \n
 ///
@@ -48,12 +51,29 @@
 /// "Topomap: A 0-dimensional homology preserving projection of high-dimensional
 /// data"\n Harish Doraiswamy, Julien Tierny, Paulo J. S. Silva, Luis Gustavo
 /// Nonato, and Claudio Silva\n Proc. of IEEE VIS 2020.\n IEEE Transactions on
-/// Visualization and Computer Graphics 27(2): 561-571, 2020.
+/// Visualization and Computer Graphics 27(2): 561-571, 2020. \n
+///
+/// "Topological Autoencoders" \n
+/// Michael Moor, Max Horn, Bastian Rieck, Karsten Borgwardt, \n
+/// Proceedings of the 37th International Conference on Machine Learning,
+/// 2020. \n
+///
+/// "Optimizing persistent homology-based functions" \n
+/// Mathieu Carriere, Frederic Chazal, Marc Glisse, Yuichi Ike,
+/// Hariprasad Kannan, Yuhei Umeda, \n
+/// Proceedings of the 38th International Conference on Machine Learning,
+/// 2021. \n
+///
+/// "Topological Autoencoders++: Fast and Accurate Cycle-Aware Dimensionality
+/// Reduction" \n
+/// Mattéo Clémot, Julie Digne, Julien Tierny, \n
+/// arXiv preprint, 2025.
 
 #pragma once
 
 #include <Debug.h>
 #include <TopoMap.h>
+#include <TopologicalDimensionReduction.h>
 
 namespace ttk {
 
@@ -78,6 +98,8 @@ namespace ttk {
       PCA = 5,
       /** TopoMap */
       TOPOMAP = 6,
+      /** AutoEncoder */
+      AE = 7,
     };
 
     inline void setSEParameters(const std::string &Affinity,
@@ -228,6 +250,9 @@ namespace ttk {
         case METHOD::TOPOMAP:
           methodName = "TopoMap (IEEE VIS 2020)";
           break;
+        case METHOD::AE:
+          methodName = "Autoencoder";
+          break;
       }
       this->printMsg("Using backend `" + methodName + "`");
     }
@@ -322,6 +347,26 @@ namespace ttk {
     bool topomap_CheckMST;
     TopoMap::STRATEGY topomap_Strategy{TopoMap::STRATEGY::KRUSKAL};
 
+    // AutoEncoder
+    bool ae_CUDA{true};
+    bool ae_Deterministic{false};
+    int ae_Seed{0};
+    int ae_Epochs{1000};
+    double ae_LearningRate{1e-2};
+    TopologicalDimensionReduction::OPTIMIZER ae_Optimizer{
+      TopologicalDimensionReduction::OPTIMIZER::ADAM};
+    TopologicalDimensionReduction::REGUL ae_Method{
+      TopologicalDimensionReduction::REGUL::ASYMMETRIC_CASCADE};
+    TopologicalDimensionReduction::MODEL ae_Model{
+      TopologicalDimensionReduction::MODEL::AUTOENCODER};
+    std::string ae_Architecture{"32 32"};
+    std::string ae_Activation{"ReLU"};
+    int ae_BatchSize{0};
+    bool ae_BatchNormalization{true};
+    double ae_RegCoefficient{1e-2};
+    bool ae_PreOptimize{false};
+    METHOD ae_PreOptimizeMethod{METHOD::PCA};
+
     // testing
     std::string ModulePath{"default"};
     std::string ModuleName{"dimensionReduction"};
@@ -334,5 +379,6 @@ namespace ttk {
     int IsDeterministic{true};
     char majorVersion_{'0'};
     bool IsInputADistanceMatrix{false};
+    bool IsInputImages{false};
   };
 } // namespace ttk

core/base/discreteMorseSandwichMPI/DiscreteMorseSandwichMPI.h

@@ -12,6 +12,7 @@
 /// "Distributed Discrete Morse Sandwich: Efficient Computation
 //  of Persistence Diagrams for Massive Scalar Data" \n
 /// Eve Le Guillou, Pierre Fortin, Julien Tierny \n
+/// IEEE Transactions on Parallel and Distributed Systems, 2025. \n
 /// arXiv:2505.21266, 2025.
 ///
 ///

core/base/persistenceDiagram/PersistenceDiagram.h

@@ -64,7 +64,8 @@
 /// "Distributed Discrete Morse Sandwich: Efficient Computation of Persistence
 /// Diagrams for Massive Scalar Data" \n
 /// Eve Le Guillou, Pierre Fortin, Julien Tierny \n
-/// https://arxiv.org/abs/2505.21266, 2025.
+/// IEEE Transactions on Parallel and Distributed Systems, 2025. \n
+/// https://arxiv.org/abs/2505.21266, 2025. \n
 /// Fast, hybrid MPI-OpenMP backend for large-scale datasets on supercomputers.
 ///
 /// \sa ttkPersistenceDiagram.cpp %for a usage example.

core/base/ripsPersistenceDiagram/RipsPersistenceDiagramUtils.h

@@ -14,7 +14,7 @@ namespace ttk::rpd {
   using value_t = double;
   constexpr value_t inf = std::numeric_limits<value_t>::infinity();
 
-  using PointCloud = std::vector<std::vector<double>>;
+  using PointCloud = std::vector<std::vector<value_t>>;
 
   using Simplex = std::vector<id_t>;
   using FiltratedSimplex = std::pair<Simplex, value_t>;

core/base/topologicalDimensionReduction/CMakeLists.txt

@@ -0,0 +1,18 @@
+ttk_add_base_library(topologicalDimensionReduction 
+  SOURCES
+    DimensionReductionModel.cpp 
+    TopologicalLoss.cpp 
+    TopologicalDimensionReduction.cpp 
+  HEADERS 
+    DimensionReductionModel.h 
+    TopologicalLoss.h 
+    TopologicalDimensionReduction.h 
+  DEPENDS
+    persistenceDiagramWarmRestartAuction 
+    ripsPersistenceDiagram)
+
+if(TTK_ENABLE_TORCH)
+target_link_libraries(topologicalDimensionReduction PUBLIC "${TORCH_LIBRARIES}")
+target_include_directories(topologicalDimensionReduction PUBLIC ${TORCH_INCLUDE_DIRS})
+target_compile_definitions(topologicalDimensionReduction PUBLIC TTK_ENABLE_TORCH) 
+endif()

core/base/topologicalDimensionReduction/DimensionReductionModel.cpp

@@ -0,0 +1,174 @@
+#include "DimensionReductionModel.h"
+#include <regex>
+
+#ifdef TTK_ENABLE_TORCH
+
+ttk::AutoEncoder::AutoEncoder(int inputDim,
+                              int latentDim,
+                              const std::string &layersDescription,
+                              const std::string &activation,
+                              bool useBN) {
+
+  std::istringstream iss(layersDescription);
+  const std::vector<std::string> hiddenDimsParsed(
+    std::istream_iterator<std::string>{iss},
+    std::istream_iterator<std::string>());
+  std::vector<unsigned> dims(1, inputDim);
+  for(const std::string &s : hiddenDimsParsed)
+    dims.push_back(std::stoi(s));
+  dims.push_back(latentDim);
+
+  const int n = dims.size() - 1;
+  encoder = torch::nn::Sequential(torch::nn::Linear(dims[0], dims[1]));
+  decoder = torch::nn::Sequential(torch::nn::Linear(dims[n], dims[n - 1]));
+  for(unsigned i = 1; i < dims.size() - 1; ++i) {
+    if(activation == "ReLU") {
+      encoder->push_back(torch::nn::ReLU());
+      decoder->push_back(torch::nn::ReLU());
+    } else if(activation == "Tanh") {
+      encoder->push_back(torch::nn::Tanh());
+      decoder->push_back(torch::nn::Tanh());
+    }
+    if(useBN) {
+      encoder->push_back(torch::nn::BatchNorm1d(dims[i]));
+      decoder->push_back(torch::nn::BatchNorm1d(dims[n - i]));
+    }
+    encoder->push_back(torch::nn::Linear(dims[i], dims[i + 1]));
+    decoder->push_back(torch::nn::Linear(dims[n - i], dims[n - (i + 1)]));
+  }
+  register_module("encoder", encoder);
+  register_module("decoder", decoder);
+}
+
+bool ttk::AutoEncoder::isStringValid(const std::string &s) {
+  return std::regex_match(s, std::regex("([0-9]+( )*)*"));
+}
+
+ttk::AutoDecoder::AutoDecoder(int inputDim,
+                              int inputSize,
+                              int latentDim,
+                              const std::string &layersDescription,
+                              const std::string &activation,
+                              bool useBN) {
+
+  std::istringstream iss(layersDescription);
+  const std::vector<std::string> hiddenDimsParsed(
+    std::istream_iterator<std::string>{iss},
+    std::istream_iterator<std::string>());
+  std::vector<unsigned> dims(1, inputDim);
+  for(const std::string &s : hiddenDimsParsed)
+    dims.push_back(std::stoi(s));
+  dims.push_back(latentDim);
+
+  latent = torch::rand({inputSize, latentDim});
+
+  const int n = dims.size() - 1;
+  decoder = torch::nn::Sequential(torch::nn::Linear(dims[n], dims[n - 1]));
+  for(unsigned i = 1; i < dims.size() - 1; ++i) {
+    if(activation == "ReLU")
+      decoder->push_back(torch::nn::ReLU());
+    else if(activation == "Tanh")
+      decoder->push_back(torch::nn::Tanh());
+    if(useBN)
+      decoder->push_back(torch::nn::BatchNorm1d(dims[n - i]));
+    decoder->push_back(torch::nn::Linear(dims[n - i], dims[n - (i + 1)]));
+  }
+  register_parameter("latent", latent, true);
+  register_module("decoder", decoder);
+}
+
+ttk::DirectOptimization::DirectOptimization(int inputSize, int latentDim) {
+  latent = torch::rand({inputSize, latentDim});
+  register_parameter("latent", latent, true);
+}
+
+ttk::ConvolutionalAutoEncoder::ConvolutionalAutoEncoder(
+  int imageSide,
+  int latentDim,
+  const std::string &layersDescription,
+  bool useBN) {
+  std::istringstream iss(layersDescription);
+  const std::vector<std::string> hiddenDimsParsed(
+    std::istream_iterator<std::string>{iss},
+    std::istream_iterator<std::string>());
+  std::vector<int> denseLayersSizes = {-1};
+  std::vector<int> convolutionalLayersChannels = {1};
+  std::vector<int> convolutionalLayersStrides;
+  for(const std::string &s : hiddenDimsParsed) {
+    if(std::regex_match(s, std::regex("c[0-9]+/[0-9]+"))) {
+      convolutionalLayersChannels.push_back(
+        std::stoi(s.substr(1, s.find('/'))));
+      convolutionalLayersStrides.push_back(
+        std::stoi(s.substr(s.find('/') + 1, s.length())));
+    } else
+      denseLayersSizes.push_back(std::stoi(s));
+  }
+  denseLayersSizes.push_back(latentDim);
+  int convolutionalOutputImageSide = imageSide;
+
+  /*** convolutional encoder ***/
+  /** unflattening for convolutional layers **/
+  encoder->push_back(torch::nn::Unflatten(
+    torch::nn::UnflattenOptions(1, {1, imageSide, imageSide})));
+  /** convolutional layers **/
+  for(unsigned c = 0; c < convolutionalLayersChannels.size() - 1; ++c) {
+    encoder->push_back(torch::nn::Conv2d(
+      torch::nn::Conv2dOptions(
+        convolutionalLayersChannels[c], convolutionalLayersChannels[c + 1], 3)
+        .padding(1)
+        .stride(convolutionalLayersStrides[c])));
+    encoder->push_back(torch::nn::ReLU());
+    convolutionalOutputImageSide /= convolutionalLayersStrides[c];
+  }
+  /** flattening for dense layers **/
+  encoder->push_back(torch::nn::Flatten());
+
+  /*** dense encoder / decoder ***/
+  denseLayersSizes[0]
+    = convolutionalLayersChannels[convolutionalLayersChannels.size() - 1]
+      * convolutionalOutputImageSide * convolutionalOutputImageSide;
+  const int n = denseLayersSizes.size() - 1;
+  encoder->push_back(
+    torch::nn::Linear(denseLayersSizes[0], denseLayersSizes[1]));
+  decoder->push_back(
+    torch::nn::Linear(denseLayersSizes[n], denseLayersSizes[n - 1]));
+  for(unsigned i = 1; i < denseLayersSizes.size() - 1; ++i) {
+    encoder->push_back(torch::nn::ReLU());
+    decoder->push_back(torch::nn::ReLU());
+    if(useBN) {
+      encoder->push_back(torch::nn::BatchNorm1d(denseLayersSizes[i]));
+      decoder->push_back(torch::nn::BatchNorm1d(denseLayersSizes[n - i]));
+    }
+    encoder->push_back(
+      torch::nn::Linear(denseLayersSizes[i], denseLayersSizes[i + 1]));
+    decoder->push_back(torch::nn::Linear(
+      denseLayersSizes[n - i], denseLayersSizes[n - (i + 1)]));
+  }
+
+  /*** convolutional decoder ***/
+  /** unflattening for convolutional layers **/
+  decoder->push_back(torch::nn::Unflatten(torch::nn::UnflattenOptions(
+    1, {convolutionalLayersChannels[convolutionalLayersChannels.size() - 1],
+        convolutionalOutputImageSide, convolutionalOutputImageSide})));
+  /** convolutional layers **/
+  for(unsigned c = convolutionalLayersChannels.size() - 1; c > 0; --c) {
+    decoder->push_back(torch::nn::ReLU());
+    decoder->push_back(torch::nn::ConvTranspose2d(
+      torch::nn::ConvTranspose2dOptions(
+        convolutionalLayersChannels[c], convolutionalLayersChannels[c - 1], 3)
+        .padding(1)
+        .stride(convolutionalLayersStrides[c - 1])
+        .output_padding(1)));
+  }
+  /** flattening for output **/
+  decoder->push_back(torch::nn::Flatten());
+
+  register_module("encoder", encoder);
+  register_module("decoder", decoder);
+}
+
+bool ttk::ConvolutionalAutoEncoder::isStringValid(const std::string &s) {
+  return std::regex_match(s, std::regex("(c[0-9]+/[0-9]+( )*)*([0-9]+( )*)*"));
+}
+
+#endif