TTK
Loading...
Searching...
No Matches
DistanceMatrixDistorsion.cpp
Go to the documentation of this file.
2#include <Geometry.h> // To check wheter a double is zero.
3#include <Os.h>
4#include <random>
5#include <vector>
6
8 this->setDebugMsgPrefix("DistanceMatrixDistorsion");
9}
10
12 const std::vector<double *> &highDistMatrix,
13 const std::vector<double *> &lowDistMatrix,
14 double &distorsionValue,
15 double *distorsionVerticesValues) const {
16 ttk::Timer timer;
17 auto n = highDistMatrix.size();
18
19 if(lowDistMatrix.size() != n) {
20 this->printErr(" Sizes mismatch: the high distance matrix has "
21 + std::to_string(n)
22 + " rows and the low distance matrix has "
23 + std::to_string(lowDistMatrix.size()) + " rows\n.");
24 return 0;
25 }
26
27 /* The computation, which is optimised for performance here, can be decomposed
28 * as follows: compute for each (x,y) delta(x,y) =
29 * (dist_low(x,y)-dist_high(x,y))^2. Compute maxi = maximum (delta(x,y)) over
30 * all (x,y). Compute delta2(x,y) = 1-delta(x,y) for each (x,y). The sim value
31 * is the mean of the n^2 values of delta2(x,y). The distorsion for a vertex
32 * x0 is the mean of delta2(x0,y) over all y's.
33 */
34
35 double maxi = 0;
36 if(distorsionVerticesValues == nullptr) {
37 this->printErr(
38 " The output pointer to the distorsionValues must be non NULL. "
39 "It must point to an allocated array of the right size.");
40 return 1;
41 }
42
43#ifdef TTK_ENABLE_OPENMP
44#pragma omp parallel for num_threads(this->threadNumber_) reduction(max \
45 : maxi) \
46 schedule(dynamic)
47#endif // TTK_ENABLE_OPENMP
48 for(size_t i = 0; i < n; i++) {
49 for(size_t j = i + 1; j < n; j++) {
50 double diff = lowDistMatrix[i][j] - highDistMatrix[i][j];
51 maxi = std::max(maxi, diff * diff);
52 }
53 }
54 const double EPS = ttk::Geometry::powIntTen(-DBL_DIG);
55 if(maxi <= EPS) // We consider maxi is equal to zero.
56 {
57 this->printMsg(
58 "The two distance matrices provided for SIM computation are equal.\n");
59 maxi = 1;
60 }
61
62 double totalSum = 0;
63
64#ifdef TTK_ENABLE_OPENMP
65#pragma omp parallel for num_threads(this->threadNumber_), reduction(+:totalSum)
66#endif // TTK_ENABLE_OPENMP
67 for(size_t i = 0; i < n; i++) {
68 double sum = 0;
69 for(size_t j = 0; j < n; j++) {
70 double diff = lowDistMatrix[i][j] - highDistMatrix[i][j];
71 double diff2 = diff * diff;
72 sum += diff2;
73 }
74 double sumHarmonized = sum / maxi;
75 distorsionVerticesValues[i] = 1 - sumHarmonized / n;
76 totalSum += 1 - sumHarmonized / n;
77 }
78
79 distorsionValue = totalSum / n;
80
81 this->printMsg("Size of output in ttk/base = " + std::to_string(n) + "\n");
82
83 this->printMsg("Computed distorsion value: "
84 + std::to_string(distorsionValue));
85 this->printMsg(ttk::debug::Separator::L2); // horizontal '-' separator
86 this->printMsg("Complete", 1, timer.getElapsedTime());
87 this->printMsg(ttk::debug::Separator::L1); // horizontal '=' separator
88
89 return 0;
90}
void setDebugMsgPrefix(const std::string &prefix)
Definition: Debug.h:364
int execute(const std::vector< double * > &highDistMatrix, const std::vector< double * > &lowDistMatrix, double &distorsionValue, double *distorsionVerticesValues) const
double getElapsedTime()
Definition: Timer.h:15
T powIntTen(const int n)
Compute the nth power of ten.
Definition: Geometry.h:360
printMsg(debug::output::GREEN+"                           "+debug::output::ENDCOLOR+debug::output::GREEN+"▒"+debug::output::ENDCOLOR+debug::output::GREEN+"▒▒▒▒▒▒▒▒▒▒▒▒▒░"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)