TTK
Loading...
Searching...
No Matches
EigenField.cpp
Go to the documentation of this file.
1#include <EigenField.h>
2#include <Laplacian.h>
3
4#if defined(TTK_ENABLE_EIGEN) && defined(TTK_ENABLE_SPECTRA)
5#include <Eigen/Eigenvalues>
6#include <Eigen/Sparse>
7
8#include <Spectra/MatOp/SparseSymMatProd.h>
9#include <Spectra/SymEigsSolver.h>
10
11#endif // TTK_ENABLE_EIGEN && TTK_ENABLE_SPECTRA
12
13// main routine
14template <typename T, class TriangulationType>
15int ttk::EigenField::execute(const TriangulationType &triangulation,
16 T *const outputFieldPointer,
17 const unsigned int eigenNumber,
18 bool computeStatistics,
19 T *const outputStatistics) const {
20
21#if defined(TTK_ENABLE_EIGEN) && defined(TTK_ENABLE_SPECTRA)
22
23 Timer tm;
24 Memory mem;
25
26 this->printMsg("Beginning computation...");
27
28#ifdef TTK_ENABLE_OPENMP
29 Eigen::setNbThreads(threadNumber_);
30#endif // TTK_ENABLE_OPENMP
31
32 using SpMat = Eigen::SparseMatrix<T>;
33 using DMat = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>;
34
35 // number of vertices
36 const auto vertexNumber = triangulation.getNumberOfVertices();
37
38 // graph laplacian of current mesh
39 SpMat lap;
40 // compute graph laplacian using cotangent weights
41 Laplacian::cotanWeights<T>(lap, *this, triangulation);
42 // lap is square
43 eigen_plain_assert(lap.cols() == lap.rows());
44
45 auto n = lap.cols();
46 auto m = eigenNumber;
47 // threshold: minimal number of eigenpairs to get a converging solution
48 const size_t minEigenNumber = 20;
49
50 if(eigenNumber == 0) {
51 // default value
52 m = n / 1000;
53 } else if(eigenNumber < minEigenNumber) {
54 m = minEigenNumber;
55 }
56
57 Spectra::SparseSymMatProd<T> op(lap);
58 Spectra::SymEigsSolver<decltype(op)> solver(op, m, 2 * m);
59
60 solver.init();
61
62 // number of eigenpairs correctly computed
63 int nconv = solver.compute(Spectra::SortRule::LargestAlge);
64
65 switch(solver.info()) {
66 case Spectra::CompInfo::NumericalIssue:
67 this->printMsg("Numerical Issue!", ttk::debug::Priority::ERROR);
68 break;
69 case Spectra::CompInfo::NotConverging:
70 this->printMsg("No Convergence! (" + std::to_string(nconv) + " out of "
71 + std::to_string(eigenNumber) + " values computed)",
73 break;
74 case Spectra::CompInfo::NotComputed:
75 this->printMsg("Invalid Input!", ttk::debug::Priority::ERROR);
76 break;
77 default:
78 break;
79 }
80
81 DMat eigenvectors = solver.eigenvectors();
82
83 auto outputEigenFunctions = static_cast<T *>(outputFieldPointer);
84
85#ifdef TTK_ENABLE_OPENMP
86#pragma omp parallel for num_threads(threadNumber_)
87#endif // TTK_ENABLE_OPENMP
88 for(SimplexId i = 0; i < vertexNumber; ++i) {
89 for(size_t j = 0; j < eigenNumber; ++j) {
90 // cannot avoid copy here...
91 outputEigenFunctions[i * eigenNumber + j] = eigenvectors(i, j);
92 }
93 }
94
95 if(computeStatistics && outputStatistics != nullptr) {
96
97 // number of statistics components
98 const int statsComp = 4;
99
100 auto outputStats = static_cast<T *>(outputStatistics);
101
102 // compute statistics on eigenfunctions
103#ifdef TTK_ENABLE_OPENMP
104#pragma omp parallel for num_threads(threadNumber_)
105#endif // TTK_ENABLE_OPENMP
106 for(SimplexId i = 0; i < vertexNumber; ++i) {
107 const auto k = i * statsComp;
108 // init current tuple computation
109 outputStats[k] = outputEigenFunctions[i * eigenNumber];
110 outputStats[k + 1] = outputEigenFunctions[i * eigenNumber];
111 outputStats[k + 2] = outputEigenFunctions[i * eigenNumber];
112 outputStats[k + 3] = outputEigenFunctions[i * eigenNumber]
113 * outputEigenFunctions[i * eigenNumber];
114 // loop from 1
115 for(size_t j = 1; j < eigenNumber; ++j) {
116 outputStats[k] = std::min<T>(
117 outputEigenFunctions[i * eigenNumber + j], outputStats[k]);
118 outputStats[k + 1] = std::max<T>(
119 outputEigenFunctions[i * eigenNumber + j], outputStats[k]);
120 outputStats[k + 2] += outputEigenFunctions[i * eigenNumber + j];
121 outputStats[k + 3] += outputEigenFunctions[i * eigenNumber + j]
122 * outputEigenFunctions[i * eigenNumber + j];
123 ;
124 }
125 }
126 }
127
128 this->printMsg("Complete", 1.0, tm.getElapsedTime(), this->threadNumber_,
129 mem.getElapsedUsage());
130
131#else
132 TTK_FORCE_USE(triangulation);
133 TTK_FORCE_USE(outputFieldPointer);
134 TTK_FORCE_USE(eigenNumber);
135 TTK_FORCE_USE(computeStatistics);
136 TTK_FORCE_USE(outputStatistics);
137
138 this->printMsg(
139 std::vector<std::string>{"Spectra support disabled, computation skipped!",
140 "Please re-compile TTK with Eigen AND Spectra "
141 "support to enable this feature."},
143#endif // TTK_ENABLE_EIGEN && TTK_ENABLE_SPECTRA
144
145 this->printMsg(ttk::debug::Separator::L1); // horizontal '=' separator
146 return 0;
147}
148
149// explicit template specializations for double and float types
150#define EIGENFIELD_SPECIALIZE(TYPE) \
151 template int ttk::EigenField::execute<TYPE>( \
152 const Triangulation &, TYPE *const, const unsigned int, bool, TYPE *const) \
153 const
154
#define TTK_FORCE_USE(x)
Force the compiler to use the function/method parameter.
Definition: BaseClass.h:57
#define EIGENFIELD_SPECIALIZE(TYPE)
Definition: EigenField.cpp:150
int threadNumber_
Definition: BaseClass.h:95
int printMsg(const std::string &msg, const debug::Priority &priority=debug::Priority::INFO, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cout) const
Definition: Debug.h:118
int execute(const TriangulationType &triangulation, T *const outputFieldPointer, const unsigned int eigenNumber=500, bool computeStatistics=false, T *const outputStatistics=nullptr) const
Definition: EigenField.cpp:15
Definition: Os.h:97
float getElapsedUsage()
Definition: Os.h:112
double getElapsedTime()
Definition: Timer.h:15
int SimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:22