TTK
Loading...
Searching...
No Matches
LowestCommonAncestor.cpp
Go to the documentation of this file.
2
3#include <algorithm>
4#include <climits>
5#include <cmath>
6#include <sstream>
7#include <stack>
8#include <string>
9
11 this->setDebugMsgPrefix("LowestCommonAncestor");
12}
13
15
16 Timer t;
17
18 // Return value
19 int retval;
20 // Eulerian transverse of the tree
21 retval = eulerianTransverse();
22 if(retval != 0) {
23 return retval;
24 }
25 // Divide the depth array in blocs and do preprocessing
26 retval = computeBlocs();
27 if(retval != 0) {
28 return retval;
29 }
30 // Preprocess the range minimum queries on the blocs
31 blocMinimumValueRMQ_.setVector(blocMinimumValue_);
32 blocMinimumValueRMQ_.setDebugLevel(debugLevel_);
33 retval = blocMinimumValueRMQ_.preprocess(true);
34 if(retval != 0) {
35 return retval;
36 }
37
38 this->printMsg("Preprocessed queries.", 1.0, t.getElapsedTime(), 1,
40
41 return 0;
42}
43
44int ttk::LowestCommonAncestor::RMQuery(const int &i, const int &j) const {
45 // Bloc of i
46 int blocI = i / blocSize_;
47 // Bloc of j
48 int blocJ = j / blocSize_;
49 // If i and j are in different blocs
50 if(blocI != blocJ) {
51 // Positions and values of the 3 range minima to compare
52 std::array<int, 3> min_pos, min_value;
53 // Position of the min in the bloc containing the ith case
54 min_pos[0] = blocI * blocSize_
55 + normalizedBlocTable_[blocToNormalizedBloc_[blocI]]
56 [i % blocSize_][blocSize_ - 1];
57 // Position of the min in the blocs between the bloc of i and j
58 min_pos[1] = ((blocJ - blocI) > 1)
59 ? blocMinimumPosition_[blocMinimumValueRMQ_.query(
60 blocI + 1, blocJ - 1)]
61 : INT_MAX;
62 // Position of the min in the bloc containing the jth case
63 min_pos[2]
64 = blocJ * blocSize_
65 + normalizedBlocTable_[blocToNormalizedBloc_[blocJ]][0][j % blocSize_];
66 // Values of the depth to compare
67 min_value[0] = nodeDepth_[min_pos[0]];
68 min_value[1] = (min_pos[1] != INT_MAX) ? nodeDepth_[min_pos[1]] : INT_MAX;
69 min_value[2] = nodeDepth_[min_pos[2]];
70 // Return the minimum of the 3
71 return min_pos[min_pos_3(min_value)];
72 } else {
73 // i and j are in the same bloc
74 return blocI * blocSize_
75 + normalizedBlocTable_[blocToNormalizedBloc_[blocI]][i % blocSize_]
76 [j % blocSize_];
77 }
78}
79
81 // Size and number of blocs
82 int sizeOfArray = static_cast<int>(nodeDepth_.size());
83 blocSize_ = static_cast<int>(log2(sizeOfArray) / 2.0);
84 int numberOfBlocs = sizeOfArray / blocSize_ + (sizeOfArray % blocSize_ != 0);
85 // Define the bloc ranges
86 for(int i = 0; i < numberOfBlocs; i++) {
87 std::pair<int, int> range;
88 range.first = i * blocSize_;
89 range.second = (i + 1) * blocSize_;
90 blocPartition_.push_back(range);
91 }
92 if(sizeOfArray % blocSize_ != 0) {
93 std::pair<int, int> range;
94 range.first = (numberOfBlocs - 1) * blocSize_;
95 range.second = sizeOfArray;
96 blocPartition_.push_back(range);
97 }
98 // Find the minimum in each bloc
99 blocMinimumValue_.resize(numberOfBlocs);
100 blocMinimumPosition_.resize(numberOfBlocs);
101 for(int i = 0; i < numberOfBlocs; i++) {
102 blocMinimumValue_[i] = nodeDepth_[blocPartition_[i].first];
103 blocMinimumPosition_[i] = blocPartition_[i].first;
104 for(int j = blocPartition_[i].first + 1; j < blocPartition_[i].second;
105 j++) {
106 if(nodeDepth_[j] < blocMinimumValue_[i]) {
107 blocMinimumValue_[i] = nodeDepth_[j];
108 blocMinimumPosition_[i] = j;
109 }
110 }
111 }
112 // Allocate the query tables
113 int numberOfTables = (1 << (blocSize_ - 1));
114 normalizedBlocTable_.resize(numberOfTables);
115 for(int i = 0; i < numberOfTables; i++) {
116 normalizedBlocTable_[i].resize(blocSize_);
117 for(int j = 0; j < blocSize_; j++) {
118 normalizedBlocTable_[i][j].resize(blocSize_);
119 }
120 }
121 // Build the query table for each possible normalized bloc
122 for(int i = 0; i < numberOfTables; i++) {
123 // Building of the ith possible bloc
124 std::vector<bool> plusOne(blocSize_ - 1);
125 int quotient = i;
126 int remain;
127 for(int j = 0; j < (blocSize_ - 1); j++) {
128 remain = quotient % 2;
129 quotient /= 2;
130 plusOne[blocSize_ - 2 - j] = static_cast<bool>(remain);
131 }
132 if(blocSize_ < 1) {
133 return -1;
134 }
135 std::vector<int> normalizedBloc(blocSize_);
136 normalizedBloc[0] = 0;
137 for(int j = 0; j < (blocSize_ - 1); j++) {
138 if(plusOne[j]) {
139 normalizedBloc[j + 1] = normalizedBloc[j] + 1;
140 } else {
141 normalizedBloc[j + 1] = normalizedBloc[j] - 1;
142 }
143 }
144 // Give the normalizedBloc to a RangeMinimumQuery object
145 RangeMinimumQuery<int> rmq(normalizedBloc);
146 rmq.setDebugLevel(debugLevel_);
147 rmq.preprocess(true);
148 // Double loop to compute all queries
149 for(int j = 0; j < blocSize_; j++) {
150 normalizedBlocTable_[i][j][j] = j;
151 for(int k = j + 1; k < blocSize_; k++) {
152 normalizedBlocTable_[i][j][k] = rmq.query(j, k);
153 normalizedBlocTable_[i][k][j] = normalizedBlocTable_[i][j][k];
154 }
155 }
156 }
157 // Determine the corresponding normalized bloc for each bloc
158 blocToNormalizedBloc_.resize(numberOfBlocs, -1);
159 for(int i = 0; i < numberOfBlocs; i++) {
160 int level = 0;
161 int tableId = 0;
162 for(int j = blocPartition_[i].first + 1; j < blocPartition_[i].second;
163 j++) {
164 if(nodeDepth_[j] > nodeDepth_[j - 1]) {
165 tableId += (1 << (blocSize_ - 2 - level));
166 }
167 level++;
168 }
169 blocToNormalizedBloc_[i] = tableId;
170 }
171 return 0;
172}
173
175 // Find the root
176 int rootId = -1;
177 for(unsigned int i = 0; i < getNumberOfNodes(); i++) {
178 if(node_[i].getAncestorId() == static_cast<int>(i)) {
179 rootId = i;
180 break;
181 }
182 }
183 if(rootId == -1) {
184 this->printErr("Tree root not found.");
185 return -1;
186 } else {
187 this->printMsg("Rout found: node id = " + std::to_string(rootId),
189 }
190 if(!(node_[rootId].getNumberOfSuccessors() > 0)) {
191 this->printErr("Tree root found with no successor.");
192 return -2;
193 }
194 // Initialize the vectors
195 nodeOrder_.clear();
196 nodeOrder_.reserve(2 * getNumberOfNodes() + 1);
197 nodeDepth_.clear();
198 nodeDepth_.reserve(2 * getNumberOfNodes() + 1);
199 nodeFirstAppearance_.clear();
200 nodeFirstAppearance_.resize(getNumberOfNodes(), -1);
201 // Transverse starting from the root
202 std::stack<int> nodeStack;
203 int depth = 0;
204 int iteration = 0;
205 std::vector<bool> isVisited(getNumberOfNodes(), false);
206 nodeStack.push(rootId);
207 while(!nodeStack.empty()) {
208 int nodeId = nodeStack.top();
209 nodeStack.pop();
210 nodeOrder_.push_back(nodeId);
211 nodeDepth_.push_back(depth);
212 if(!isVisited[nodeId]) {
213 // Add ancestor
214 if(nodeId != rootId) {
215 nodeStack.push(node_[nodeId].getAncestorId());
216 }
217 // Add successors
218 int numberOfSuccessors = node_[nodeId].getNumberOfSuccessors();
219 for(int i = 0; i < numberOfSuccessors; i++) {
220 nodeStack.push(node_[nodeId].getSuccessorId(i));
221 }
222 nodeFirstAppearance_[nodeId] = iteration;
223 isVisited[nodeId] = true;
224 }
225 // Next depth
226 if(!nodeStack.empty()) {
227 if(nodeStack.top() == node_[nodeId].getAncestorId()) {
228 depth--;
229 } else {
230 depth++;
231 }
232 }
233 iteration++;
234 }
235 return 0;
236}
void setDebugMsgPrefix(const std::string &prefix)
Definition: Debug.h:364
virtual int setDebugLevel(const int &debugLevel)
Definition: Debug.cpp:147
int RMQuery(const int &i, const int &j) const
Class to answer range minimum queries in an array in constant time after a linearithmic time preproce...
int preprocess(const bool silent=false)
int query(int i, int j) const
double getElapsedTime()
Definition: Timer.h:15
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)