TTK
Loading...
Searching...
No Matches
RipsComplex.cpp
Go to the documentation of this file.
1#include <RipsComplex.h>
2
3#include <limits>
4
6 this->setDebugMsgPrefix("RipsComplex");
7}
8
9template <size_t n>
10struct LocCell {
11 double diam;
12 std::array<ttk::SimplexId, n> verts;
13};
14
15static void computeEdges(std::vector<ttk::SimplexId> &connectivity,
16 std::vector<double> &diameters,
17 const double epsilon,
18 const std::vector<std::vector<double>> &distanceMatrix,
19 const int nThreads) {
20
21 TTK_FORCE_USE(nThreads);
22
23 std::vector<std::vector<LocCell<1>>> edges(distanceMatrix.size());
24
25#ifdef TTK_ENABLE_OPENMP
26#pragma omp parallel for num_threads(nThreads) schedule(dynamic)
27#endif // TTK_ENABLE_OPENMP
28 for(size_t i = 0; i < distanceMatrix.size(); ++i) {
29 for(size_t j = i + 1; j < distanceMatrix.size(); ++j) {
30 if(distanceMatrix[i][j] < epsilon) {
31 edges[i].emplace_back(LocCell<1>{
32 distanceMatrix[i][j],
33 std::array<ttk::SimplexId, 1>{static_cast<ttk::SimplexId>(j)}});
34 }
35 }
36 }
37
38 std::vector<size_t> psum(edges.size() + 1);
39 for(size_t i = 0; i < edges.size(); ++i) {
40 psum[i + 1] = psum[i] + edges[i].size();
41 }
42
43 const auto nCells{psum.back()};
44 diameters.resize(nCells);
45 connectivity.resize(2 * nCells);
46
47#ifdef TTK_ENABLE_OPENMP
48#pragma omp parallel for num_threads(nThreads) schedule(dynamic)
49#endif // TTK_ENABLE_OPENMP
50 for(size_t i = 0; i < edges.size(); ++i) {
51 for(size_t j = 0; j < edges[i].size(); ++j) {
52 connectivity[2 * (psum[i] + j) + 0] = static_cast<ttk::SimplexId>(i);
53 connectivity[2 * (psum[i] + j) + 1] = edges[i][j].verts[0];
54 diameters[psum[i] + j] = edges[i][j].diam;
55 }
56 }
57}
58
59static inline void maxAssign(double &a, const double b) {
60 if(a < b) {
61 a = b;
62 }
63}
64
65static void
66 computeTriangles(std::vector<ttk::SimplexId> &connectivity,
67 std::vector<double> &diameters,
68 const double epsilon,
69 const std::vector<std::vector<double>> &distanceMatrix,
70 const int nThreads) {
71
72 TTK_FORCE_USE(nThreads);
73
74 std::vector<std::vector<LocCell<2>>> triangles(distanceMatrix.size());
75
76#ifdef TTK_ENABLE_OPENMP
77#pragma omp parallel for num_threads(nThreads) schedule(dynamic)
78#endif // TTK_ENABLE_OPENMP
79 for(size_t i = 0; i < distanceMatrix.size(); ++i) {
80 for(size_t j = i + 1; j < distanceMatrix.size(); ++j) {
81 if(distanceMatrix[i][j] > epsilon) {
82 continue;
83 }
84 for(size_t k = j + 1; k < distanceMatrix.size(); ++k) {
85 if(distanceMatrix[i][k] > epsilon || distanceMatrix[j][k] > epsilon) {
86 continue;
87 }
88 auto diam{distanceMatrix[i][j]};
89 maxAssign(diam, distanceMatrix[i][k]);
90 maxAssign(diam, distanceMatrix[j][k]);
91 triangles[i].emplace_back(
92 LocCell<2>{diam, std::array<ttk::SimplexId, 2>{
93 static_cast<ttk::SimplexId>(j),
94 static_cast<ttk::SimplexId>(k),
95 }});
96 }
97 }
98 }
99
100 std::vector<size_t> psum(triangles.size() + 1);
101 for(size_t i = 0; i < triangles.size(); ++i) {
102 psum[i + 1] = psum[i] + triangles[i].size();
103 }
104
105 const auto nCells{psum.back()};
106 diameters.resize(nCells);
107 connectivity.resize(3 * nCells);
108
109#ifdef TTK_ENABLE_OPENMP
110#pragma omp parallel for num_threads(nThreads) schedule(dynamic)
111#endif // TTK_ENABLE_OPENMP
112 for(size_t i = 0; i < triangles.size(); ++i) {
113 for(size_t j = 0; j < triangles[i].size(); ++j) {
114 connectivity[3 * (psum[i] + j) + 0] = static_cast<ttk::SimplexId>(i);
115 connectivity[3 * (psum[i] + j) + 1] = triangles[i][j].verts[0];
116 connectivity[3 * (psum[i] + j) + 2] = triangles[i][j].verts[1];
117 diameters[psum[i] + j] = triangles[i][j].diam;
118 }
119 }
120}
121
122static void
123 computeTetras(std::vector<ttk::SimplexId> &connectivity,
124 std::vector<double> &diameters,
125 const double epsilon,
126 const std::vector<std::vector<double>> &distanceMatrix,
127 const int nThreads) {
128
129 TTK_FORCE_USE(nThreads);
130
131 std::vector<std::vector<LocCell<3>>> tetras(distanceMatrix.size());
132
133#ifdef TTK_ENABLE_OPENMP
134#pragma omp parallel for num_threads(nThreads) schedule(dynamic)
135#endif // TTK_ENABLE_OPENMP
136 for(size_t i = 0; i < distanceMatrix.size(); ++i) {
137 for(size_t j = i + 1; j < distanceMatrix.size(); ++j) {
138 if(distanceMatrix[i][j] > epsilon) {
139 continue;
140 }
141 for(size_t k = j + 1; k < distanceMatrix.size(); ++k) {
142 if(distanceMatrix[i][k] > epsilon || distanceMatrix[j][k] > epsilon) {
143 continue;
144 }
145 for(size_t l = k + 1; l < distanceMatrix.size(); ++l) {
146 if(distanceMatrix[i][l] > epsilon || distanceMatrix[j][l] > epsilon
147 || distanceMatrix[k][l] > epsilon) {
148 continue;
149 }
150 auto diam{distanceMatrix[i][j]};
151 maxAssign(diam, distanceMatrix[i][k]);
152 maxAssign(diam, distanceMatrix[i][l]);
153 maxAssign(diam, distanceMatrix[j][k]);
154 maxAssign(diam, distanceMatrix[j][l]);
155 maxAssign(diam, distanceMatrix[k][l]);
156 tetras[i].emplace_back(
157 LocCell<3>{diam, std::array<ttk::SimplexId, 3>{
158 static_cast<ttk::SimplexId>(j),
159 static_cast<ttk::SimplexId>(k),
160 static_cast<ttk::SimplexId>(l),
161 }});
162 }
163 }
164 }
165 }
166
167 std::vector<size_t> psum(tetras.size() + 1);
168 for(size_t i = 0; i < tetras.size(); ++i) {
169 psum[i + 1] = psum[i] + tetras[i].size();
170 }
171
172 const auto nCells{psum.back()};
173 diameters.resize(nCells);
174 connectivity.resize(4 * nCells);
175
176#ifdef TTK_ENABLE_OPENMP
177#pragma omp parallel for num_threads(nThreads) schedule(dynamic)
178#endif // TTK_ENABLE_OPENMP
179 for(size_t i = 0; i < tetras.size(); ++i) {
180 for(size_t j = 0; j < tetras[i].size(); ++j) {
181 connectivity[4 * (psum[i] + j) + 0] = static_cast<ttk::SimplexId>(i);
182 connectivity[4 * (psum[i] + j) + 1] = tetras[i][j].verts[0];
183 connectivity[4 * (psum[i] + j) + 2] = tetras[i][j].verts[1];
184 connectivity[4 * (psum[i] + j) + 3] = tetras[i][j].verts[2];
185 diameters[psum[i] + j] = tetras[i][j].diam;
186 }
187 }
188}
189
191 double *const density,
192 const std::vector<std::vector<double>> &distanceMatrix) const {
193
194 const auto sq = [](const double a) -> double { return a * a; };
195
196#ifdef TTK_ENABLE_OPENMP
197#pragma omp parallel for num_threads(this->threadNumber_)
198#endif // TTK_ENABLE_OPENMP
199 for(size_t i = 0; i < distanceMatrix.size(); ++i) {
200 density[i] = 0.0;
201 for(size_t j = 0; j < distanceMatrix.size(); ++j) {
202 if(i == j) {
203 density[i] += 1.0;
204 }
205 density[i]
206 += std::exp(-sq(distanceMatrix[i][j]) / (2.0 * sq(this->StdDev)));
207 }
208 }
209
210 return 0;
211}
212
214 const SimplexId nPoints,
215 std::array<double *const, 3> diamStats,
216 const std::vector<SimplexId> &connectivity,
217 const std::vector<double> &cellDiameters) const {
218
219 const auto nCells{cellDiameters.size()};
220 if(nCells != connectivity.size() / (this->OutputDimension + 1)) {
221 this->printErr("Cell number mismatch");
222 return -1;
223 }
224
225 std::vector<size_t> nCellsAroundVert(nPoints, 0);
226
227#ifdef TTK_ENABLE_OPENMP
228#pragma omp parallel for num_threads(this->threadNumber_)
229#endif // TTK_ENABLE_OPENMP
230 for(SimplexId i = 0; i < nPoints; ++i) {
231 diamStats[0][i] = this->Epsilon; // min
232 diamStats[1][i] = 0.0; // mean
233 diamStats[2][i] = 0.0; // max
234 }
235
236 for(size_t i = 0; i < nCells; ++i) {
237 for(int j = 0; j < this->OutputDimension + 1; ++j) {
238 const auto p{connectivity[i * (this->OutputDimension + 1) + j]};
239 nCellsAroundVert[p]++;
240 diamStats[0][p] = std::min(diamStats[0][p], cellDiameters[i]);
241 diamStats[1][p] += cellDiameters[i];
242 diamStats[2][p] = std::max(diamStats[2][p], cellDiameters[i]);
243 }
244 }
245
246#ifdef TTK_ENABLE_OPENMP
247#pragma omp parallel for num_threads(this->threadNumber_)
248#endif // TTK_ENABLE_OPENMP
249 for(SimplexId i = 0; i < nPoints; ++i) {
250 if(nCellsAroundVert[i] == 0) {
251 diamStats[0][i] = 0.0; // min
252 diamStats[1][i] = 0.0; // mean
253 diamStats[2][i] = 0.0; // max
254 } else {
255 diamStats[1][i] /= nCellsAroundVert[i];
256 }
257 }
258
259 return 0;
260}
261
263 std::vector<SimplexId> &connectivity,
264 std::vector<double> &diameters,
265 std::array<double *const, 3> diamStats,
266 const std::vector<std::vector<double>> &distanceMatrix,
267 double *const density) const {
268
269 Timer tm{};
270
271 if(distanceMatrix.empty()
272 || distanceMatrix.size() != distanceMatrix[0].size()) {
273 this->printErr("Invalid distance matrix");
274 return 1;
275 }
276
277 Timer tm_rips{};
278
279 if(this->OutputDimension == 1) {
280 computeEdges(connectivity, diameters, this->Epsilon, distanceMatrix,
281 this->threadNumber_);
282 } else if(this->OutputDimension == 2) {
283 computeTriangles(connectivity, diameters, this->Epsilon, distanceMatrix,
284 this->threadNumber_);
285 } else if(this->OutputDimension == 3) {
286 computeTetras(connectivity, diameters, this->Epsilon, distanceMatrix,
287 this->threadNumber_);
288 }
289
290 this->printMsg("Generated Rips complex from distance matrix", 1.0,
291 tm_rips.getElapsedTime(), this->threadNumber_,
293
294 this->computeDiameterStats(
295 distanceMatrix.size(), diamStats, connectivity, diameters);
296
297 if(this->ComputeGaussianDensity) {
298 this->computeGaussianDensity(density, distanceMatrix);
299 }
300
301 this->printMsg("Complete", 1.0, tm.getElapsedTime(), this->threadNumber_);
302 return 0;
303}
#define TTK_FORCE_USE(x)
Force the compiler to use the function/method parameter.
Definition BaseClass.h:57
void setDebugMsgPrefix(const std::string &prefix)
Definition Debug.h:364
int computeGaussianDensity(double *const density, const std::vector< std::vector< double > > &distanceMatrix) const
Compute Gaussian density on points.
int computeDiameterStats(const SimplexId nPoints, std::array< double *const, 3 > diamStats, const std::vector< SimplexId > &connectivity, const std::vector< double > &cellDiameters) const
Compute diameter statistics on points.
int execute(std::vector< SimplexId > &connectivity, std::vector< double > &diameters, std::array< double *const, 3 > diamStats, const std::vector< std::vector< double > > &distanceMatrix, double *const density=nullptr) const
Main entry point.
int SimplexId
Identifier type for simplices of any dimension.
Definition DataTypes.h:22
double diam
std::array< ttk::SimplexId, n > verts
printMsg(debug::output::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/|__ _|"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)