1 | #include "mpi.hpp" |
---|
2 | #include <cstdlib> |
---|
3 | #include <cmath> |
---|
4 | #include <iostream> |
---|
5 | #include <cassert> |
---|
6 | #include "tree.hpp" |
---|
7 | #include "elt.hpp" |
---|
8 | #include "intersect.hpp" |
---|
9 | #include <vector> |
---|
10 | #include <set> |
---|
11 | #include <algorithm> |
---|
12 | |
---|
13 | #include "node.hpp" |
---|
14 | |
---|
15 | namespace sphereRemap { |
---|
16 | |
---|
17 | |
---|
18 | #define UPDATE_EVERY 1 |
---|
19 | |
---|
20 | using namespace std; |
---|
21 | |
---|
22 | //#ifdef DEBUG |
---|
23 | //std::map<void*, struct mem_info> _mem_deb; // reference counter |
---|
24 | // |
---|
25 | //void NodePtr::unlink() |
---|
26 | //{ |
---|
27 | // _mem_deb[ptr].ref_cnt -= 1; |
---|
28 | // if (_mem_deb[ptr].ref_cnt == 0) // we were last reference |
---|
29 | // { |
---|
30 | // if (_mem_deb[ptr].stat == DELETED) // delete has already been called, everything is fine, free memory now |
---|
31 | // { |
---|
32 | // _mem_deb.erase(ptr); |
---|
33 | // ::delete [] ((char *) ptr); |
---|
34 | // } |
---|
35 | // else // no more pointer, but memory not freed -> memory leak !! |
---|
36 | // { |
---|
37 | // std::cerr << "WARNING: Memory leak created at address " << ptr << "."; |
---|
38 | // assert(_mem_deb[ptr].stat == ALLOCATED); // and not BORROWED |
---|
39 | //#ifdef AUTO_LEAK_FIX |
---|
40 | // // free automatically, just for fun |
---|
41 | // _mem_deb.erase(ptr); |
---|
42 | // ptr->~Node(); // call destructor, since `Node::delete` has not been applied it has not been called |
---|
43 | // ::delete [] ((char *) ptr); |
---|
44 | // std::cerr << " LEAK FIXED."; |
---|
45 | //#endif |
---|
46 | // |
---|
47 | // std::cerr << std::endl; |
---|
48 | // } |
---|
49 | // } |
---|
50 | //} |
---|
51 | // |
---|
52 | //void memory_report() |
---|
53 | //{ |
---|
54 | // for (std::map<void*, struct mem_info>::iterator it = _mem_deb.begin(); it != _mem_deb.end(); it++) |
---|
55 | // { |
---|
56 | // if (it->second.stat == BORROWED) continue; |
---|
57 | // std::cerr << "Node at " << it->first << " has " << it->second.ref_cnt << " references to it and " |
---|
58 | // << ((it->second.stat == DELETED) ? "has" : "has *not*") << " been deleted." << std::endl; |
---|
59 | //#ifdef AUTO_LEAK_FIX |
---|
60 | // if (it->second.stat == ALLOCATED) // `Node::delete` has not been called |
---|
61 | // ((Node *) it->first)->~Node(); |
---|
62 | // delete [] ((char *) it->first); |
---|
63 | //#endif |
---|
64 | // } |
---|
65 | //} |
---|
66 | // |
---|
67 | //int ref_cnt(void *ptr) |
---|
68 | //{ |
---|
69 | // return _mem_deb[ptr].ref_cnt; |
---|
70 | //} |
---|
71 | //#endif |
---|
72 | |
---|
73 | int compareDist(NodePtr n1, NodePtr n2) |
---|
74 | { |
---|
75 | return squaredist(n1->tree->ref->centre, n1->centre) |
---|
76 | < squaredist(n2->tree->ref->centre, n2->centre); |
---|
77 | } |
---|
78 | |
---|
79 | /* On level `level` find the node in our subtree that is closest to `src` and return through argument `closest`. |
---|
80 | The return value is just for recursive calling */ |
---|
81 | void Node::findClosest(int level, NodePtr src, double& minDist2, NodePtr &closest) |
---|
82 | { |
---|
83 | double r2; |
---|
84 | |
---|
85 | r2 = squaredist(src->centre, this->centre); |
---|
86 | if (level == this->level) |
---|
87 | { |
---|
88 | if (r2 < minDist2 || closest == NULL) |
---|
89 | { |
---|
90 | minDist2 = r2; |
---|
91 | closest = this; |
---|
92 | } |
---|
93 | } |
---|
94 | else if (r2 < radius*radius) |
---|
95 | { |
---|
96 | for(int i = 0; i < child.size(); i++) |
---|
97 | child[i]->findClosest(level, src, minDist2, closest); |
---|
98 | } |
---|
99 | } |
---|
100 | |
---|
101 | NodePtr Node::farthest(vector<NodePtr >& list) |
---|
102 | { |
---|
103 | assert(this); |
---|
104 | double dmax = -HUGE_VAL; |
---|
105 | NodePtr node = NULL; |
---|
106 | for (size_t i=0; i < list.size(); i++) |
---|
107 | { |
---|
108 | double d = ds(this->centre, list[i]->centre); |
---|
109 | if (d > dmax) { |
---|
110 | node = list[i]; |
---|
111 | dmax = d; |
---|
112 | } |
---|
113 | } |
---|
114 | return node; |
---|
115 | } |
---|
116 | |
---|
117 | /* returns element in `list` cosest to us (`this`) by meassure `ds`. |
---|
118 | If `n` is negative finds element furthest away instead. */ |
---|
119 | NodePtr Node::closest(vector<NodePtr >& list, int n) |
---|
120 | { |
---|
121 | assert(this); |
---|
122 | double dmin2 = (n>0) ? HUGE_VAL : -HUGE_VAL; |
---|
123 | NodePtr node = NULL; |
---|
124 | for (int i = 0; i < list.size(); i++) |
---|
125 | { |
---|
126 | double d2 = squaredist(this->centre, list[i]->centre); |
---|
127 | if (n * (d2 - dmin2) < 0) |
---|
128 | { |
---|
129 | node = list[i]; |
---|
130 | dmin2 = d2; |
---|
131 | } |
---|
132 | } |
---|
133 | return node; |
---|
134 | } |
---|
135 | |
---|
136 | // make sure we will be able to accomodate `node` |
---|
137 | void Node::move(const NodePtr node) // this->leafCount may be 0 |
---|
138 | { |
---|
139 | double w = double(node->leafCount)/double(node->leafCount + this->leafCount); |
---|
140 | Coord oldCentre = this->centre; |
---|
141 | this->centre = proj(this->centre * (1-w) + node->centre * w); |
---|
142 | this->leafCount += node->leafCount; |
---|
143 | this->radius += arcdist(oldCentre, this->centre) + 1e-9; |
---|
144 | } |
---|
145 | |
---|
146 | void Node::remove(const NodePtr node) |
---|
147 | { |
---|
148 | if (&node == NULL) return; |
---|
149 | double w = double(node->leafCount) / this->leafCount; |
---|
150 | Coord newCentre = proj(this->centre - node->centre * w); |
---|
151 | this->leafCount -= node->leafCount; |
---|
152 | this->radius += arcdist(newCentre, this->centre) + 1e-9; |
---|
153 | this->centre = newCentre; |
---|
154 | this->updateCount++; |
---|
155 | } |
---|
156 | |
---|
157 | void Node::inflate(const NodePtr node) |
---|
158 | { |
---|
159 | double d = arcdist(this->centre, node->centre); |
---|
160 | double r = node->radius; |
---|
161 | if (this->radius < d + r) |
---|
162 | this->radius = d + r + 1e-9; |
---|
163 | } |
---|
164 | |
---|
165 | /* recomputes the radius which currently might be much larger than necessary */ |
---|
166 | void Node::update() |
---|
167 | { |
---|
168 | Coord centre = ORIGIN; |
---|
169 | int n = 0; |
---|
170 | for (int i = 0; i < this->child.size(); i++) { |
---|
171 | centre = centre + this->child[i]->centre * this->child[i]->leafCount; |
---|
172 | n += this->child[i]->leafCount; |
---|
173 | } |
---|
174 | this->centre = proj(centre); |
---|
175 | this->leafCount = n; |
---|
176 | |
---|
177 | double R = 0; |
---|
178 | for (int i = 0; i < this->child.size(); i++) |
---|
179 | { |
---|
180 | double d = arcdist(this->centre, this->child[i]->centre); |
---|
181 | double r = this->child[i]->radius; |
---|
182 | if (R < d + r) R = d + r; |
---|
183 | } |
---|
184 | this->radius = R + 1e-9; |
---|
185 | this->updateCount = 0; |
---|
186 | if (child.size()) |
---|
187 | level = child[0]->level + 1; |
---|
188 | } |
---|
189 | |
---|
190 | void Node::output(ostream& flux, int level, int color) |
---|
191 | { |
---|
192 | if (level==this->level) |
---|
193 | { |
---|
194 | flux<<centre.x<<" , "<<centre.y<<" , "<<centre.z<<" , "<<radius<<endl ; |
---|
195 | } |
---|
196 | else |
---|
197 | { |
---|
198 | for (int i = 0; i < this->child.size(); i++) child[i]->output(flux,level,color) ; |
---|
199 | } |
---|
200 | } |
---|
201 | |
---|
202 | |
---|
203 | /* void Node::append(NodePtr node); |
---|
204 | { |
---|
205 | if (node->level == this->level -1) |
---|
206 | { |
---|
207 | // new node is one level lower (correct level) |
---|
208 | this->child.append(node); |
---|
209 | node->parent = this; |
---|
210 | } |
---|
211 | else if (node->level < this->level -1) |
---|
212 | { |
---|
213 | // if new node is on even lower level, recursively append to closest child |
---|
214 | node->closest(this->child)->append(node); |
---|
215 | } |
---|
216 | else |
---|
217 | { |
---|
218 | cerr << "Error: Attempted node insertion with invalid level." << endl; |
---|
219 | exit(1); |
---|
220 | } |
---|
221 | }*/ |
---|
222 | |
---|
223 | bool find_in_tree1(Node* node) |
---|
224 | { |
---|
225 | if (node == node->tree->root) return true; |
---|
226 | if (node->parent == NULL) |
---|
227 | { |
---|
228 | cerr << "Cannot find!" << endl; |
---|
229 | return false; |
---|
230 | } |
---|
231 | return find_in_tree1(node->parent); |
---|
232 | } |
---|
233 | |
---|
234 | bool find_in_tree2(NodePtr node, NodePtr ref) |
---|
235 | { |
---|
236 | for (int i = 0; i < ref->child.size(); i++) |
---|
237 | { |
---|
238 | if (node == ref->child[i]) |
---|
239 | { |
---|
240 | cerr << "find2: " << ref << " -> " << ref->child[i] << endl; |
---|
241 | return true; |
---|
242 | } |
---|
243 | else if (find_in_tree2(node, ref->child[i])) |
---|
244 | { |
---|
245 | cerr << "find2: " << ref << " -> " << ref->child[i] << endl; |
---|
246 | return true; |
---|
247 | } |
---|
248 | } |
---|
249 | return false; |
---|
250 | } |
---|
251 | |
---|
252 | |
---|
253 | /* This appends `this` to the node `node` */ |
---|
254 | NodePtr insert(NodePtr thIs, NodePtr node) |
---|
255 | { |
---|
256 | int la = thIs->level; // node to be inserted |
---|
257 | int lb = node->level; // node where insertation |
---|
258 | assert(la < lb); // node to be inserted must have lower level then parent |
---|
259 | //if (thIs->parent) assert(find_in_tree1(thIs) == true); |
---|
260 | NodePtr q = NULL; |
---|
261 | NodePtr chd = NULL; |
---|
262 | node->move(thIs); |
---|
263 | if (la == lb - 1) |
---|
264 | { |
---|
265 | node->child.push_back(thIs); |
---|
266 | thIs->parent = node; |
---|
267 | if (node->child.size() > MAX_NODE_SZ && node->tree->canSplit() ) // with us as additional child `node` is now too large :( |
---|
268 | return (node->reinserted || node->parent == NULL) ? split(node) : reinsert(node); |
---|
269 | } |
---|
270 | else // la < lb - 1 |
---|
271 | { |
---|
272 | chd = thIs->closest(node->child); |
---|
273 | q = insert(thIs, chd); |
---|
274 | } |
---|
275 | if ((node->updateCount + 1) % UPDATE_EVERY == 0) |
---|
276 | node->update(); |
---|
277 | else |
---|
278 | { |
---|
279 | if (q) node->remove(q); |
---|
280 | node->inflate(chd); |
---|
281 | } |
---|
282 | |
---|
283 | |
---|
284 | return q; |
---|
285 | } |
---|
286 | |
---|
287 | typedef NodePtr pNode; |
---|
288 | |
---|
289 | /* move `frac` of our children to a new node which is returned */ |
---|
290 | NodePtr reinsert(NodePtr thIs) |
---|
291 | { |
---|
292 | thIs->tree->ref = thIs; |
---|
293 | std::sort(thIs->child.begin(), thIs->child.end(), compareDist); |
---|
294 | int out = (int) (frac*thIs->child.size()); |
---|
295 | |
---|
296 | /* make sure out is only so big that there are still MIN_NODE_SZ children after removing out */ |
---|
297 | if (thIs->child.size() - out < MIN_NODE_SZ) out = thIs->child.size() - MIN_NODE_SZ; |
---|
298 | |
---|
299 | |
---|
300 | /* transfere out children from us to a new node q which will be returned */ |
---|
301 | NodePtr q = new Node; |
---|
302 | q->tree = thIs->tree; |
---|
303 | q->child.resize(out); |
---|
304 | for (int i = thIs->child.size() - out; i < thIs->child.size(); i++) |
---|
305 | { |
---|
306 | thIs->tree->push_back(thIs->child[i]); |
---|
307 | int k = i - (thIs->child.size() - out); |
---|
308 | q->child[k] = thIs->child[i]; |
---|
309 | q->child[k]->parent = q; |
---|
310 | } |
---|
311 | thIs->child.resize(thIs->child.size() - out); |
---|
312 | thIs->update(); |
---|
313 | q->update(); |
---|
314 | thIs->reinserted = true; // avoid infinite loop of reinserting the same node, by marking it as reinserted and stop if same node arrives at same place again |
---|
315 | thIs->tree->ri = 1; |
---|
316 | |
---|
317 | return q; |
---|
318 | } |
---|
319 | |
---|
320 | /* move around nodes that are far away from the centre of their parents in order reduce radia of the circles |
---|
321 | leading to faster look-up times because of less redundancies between nodes. |
---|
322 | TODO cite paper for Slim SS-tree */ |
---|
323 | void slim2(NodePtr thIs, int level, int minNodeSize) |
---|
324 | { |
---|
325 | bool out; |
---|
326 | double distChild; |
---|
327 | |
---|
328 | #ifdef DEBUG |
---|
329 | // assert(ref_cnt(thIs) >= thIs->child.size() + 1 /*parent*/ + 1 /*thIs*/); |
---|
330 | #endif |
---|
331 | if (thIs->level==level) |
---|
332 | { |
---|
333 | /* |
---|
334 | out = false; |
---|
335 | while (!out) |
---|
336 | { |
---|
337 | // remove child which is farthest away from the centre and try to reinsert it into the tree |
---|
338 | double distMax = 0; |
---|
339 | int itMax = -1; |
---|
340 | |
---|
341 | for (int i = 0; i < thIs->child.size(); i++) |
---|
342 | { |
---|
343 | distChild = arcdist(thIs->centre, thIs->child[i]->centre) + thIs->child[i]->radius; |
---|
344 | if (distChild > distMax) |
---|
345 | { |
---|
346 | distMax = distChild; |
---|
347 | itMax = i; |
---|
348 | } |
---|
349 | } |
---|
350 | if (transferNode(thIs->tree->root, thIs, thIs->child[itMax])) |
---|
351 | { |
---|
352 | thIs->child.erase(thIs->child.begin()+itMax); |
---|
353 | out = false; |
---|
354 | } |
---|
355 | else |
---|
356 | out = true; |
---|
357 | |
---|
358 | if (thIs->child.size() < minNodeSize) out = true; |
---|
359 | } |
---|
360 | */ |
---|
361 | if (thIs->tree-> isActiveOkSplit && thIs->tree->levelSize[thIs->tree->assignLevel] <= thIs->tree->keepNodes) |
---|
362 | { |
---|
363 | |
---|
364 | return ; |
---|
365 | } |
---|
366 | for (int i = 0; i < thIs->child.size(); i++) |
---|
367 | { |
---|
368 | std::vector<NodePtr> before; |
---|
369 | if (transferNode(thIs->tree->root, thIs, thIs->child[i])) |
---|
370 | { |
---|
371 | before=thIs->child ; |
---|
372 | thIs->child.erase(thIs->child.begin()+i); |
---|
373 | i--; |
---|
374 | } |
---|
375 | } |
---|
376 | |
---|
377 | |
---|
378 | if (thIs->child.size() < minNodeSize && thIs->level < thIs->tree->root->level) |
---|
379 | { |
---|
380 | thIs->tree->decreaseLevelSize(thIs->level); |
---|
381 | for(int i = 0; i < thIs->child.size(); i++) |
---|
382 | thIs->tree->push_back(thIs->child[i]); |
---|
383 | thIs->child.resize(0); |
---|
384 | } |
---|
385 | else thIs->update(); |
---|
386 | } |
---|
387 | else |
---|
388 | { |
---|
389 | int newChildCount = 0; |
---|
390 | for (int i = 0; i < thIs->child.size(); i++) |
---|
391 | { |
---|
392 | if (thIs == thIs->tree->root) |
---|
393 | { |
---|
394 | // keep at least one child for root |
---|
395 | if (i == thIs->child.size()-1 && newChildCount == 0) |
---|
396 | { |
---|
397 | thIs->child[newChildCount] = thIs->child[i]; |
---|
398 | newChildCount++; |
---|
399 | break; |
---|
400 | } |
---|
401 | } |
---|
402 | |
---|
403 | slim2(thIs->child[i], level); |
---|
404 | if (thIs->child[i]->child.size() != 0) // thIs->child[i] is not a leaf |
---|
405 | { |
---|
406 | thIs->child[newChildCount] = thIs->child[i]; |
---|
407 | newChildCount++; |
---|
408 | } /* else FIXME sometimes this child must be deleted (otherwise leak) sometimes not (otherwise segfault) |
---|
409 | maybe delete not here but when transfered |
---|
410 | delete thIs->child[i]; // if our child does not make any grand-children what good is it? -> remove! |
---|
411 | */ |
---|
412 | |
---|
413 | } |
---|
414 | thIs->child.resize(newChildCount); |
---|
415 | |
---|
416 | if (thIs->child.size() < MIN_NODE_SZ && thIs->level < thIs->tree->root->level) |
---|
417 | { |
---|
418 | thIs->tree->decreaseLevelSize(thIs->level); |
---|
419 | for (int i = 0; i < thIs->child.size(); i++) |
---|
420 | thIs->tree->push_front(thIs->child[i]); |
---|
421 | thIs->child.resize(0); |
---|
422 | } |
---|
423 | else thIs->update(); |
---|
424 | } |
---|
425 | |
---|
426 | } |
---|
427 | |
---|
428 | bool transferNode(NodePtr thIs, NodePtr parent, NodePtr node) |
---|
429 | { |
---|
430 | if (parent == thIs) return false; |
---|
431 | |
---|
432 | if (thIs->level == parent->level) |
---|
433 | { |
---|
434 | if ( (thIs->child.size() < MAX_NODE_SZ || thIs->tree->isActiveOkSplit) && thIs->child.size() >= MIN_NODE_SZ) |
---|
435 | { |
---|
436 | insert(node, thIs); |
---|
437 | return true; |
---|
438 | } |
---|
439 | else |
---|
440 | return false; |
---|
441 | } |
---|
442 | else |
---|
443 | { |
---|
444 | for (int i = 0; i < thIs->child.size(); i++) |
---|
445 | { |
---|
446 | if (arcdist(thIs->child[i]->centre, node->centre) + node->radius < thIs->child[i]->radius) |
---|
447 | { |
---|
448 | if (transferNode(thIs->child[i], parent, node)) |
---|
449 | { |
---|
450 | thIs->update(); |
---|
451 | return true; |
---|
452 | } |
---|
453 | } |
---|
454 | } |
---|
455 | return false; |
---|
456 | } |
---|
457 | } |
---|
458 | |
---|
459 | |
---|
460 | |
---|
461 | NodePtr split(NodePtr thIs) |
---|
462 | { |
---|
463 | thIs->tree->increaseLevelSize(thIs->level); |
---|
464 | NodePtr p = new Node; |
---|
465 | NodePtr q = new Node; |
---|
466 | p->level = q->level = thIs->level; |
---|
467 | p->reinserted = q->reinserted = false; |
---|
468 | p->parent = q->parent = thIs->parent; |
---|
469 | p->tree = q->tree = thIs->tree; |
---|
470 | |
---|
471 | |
---|
472 | p->child.resize(MAX_NODE_SZ/2); |
---|
473 | q->child.resize(MAX_NODE_SZ/2 + 1); |
---|
474 | assert(thIs->child.size() == MAX_NODE_SZ+1); |
---|
475 | thIs->tree->ref = thIs->closest(thIs->child, FARTHEST); // farthest from centre |
---|
476 | std::sort(thIs->child.begin(), thIs->child.end(), compareDist); |
---|
477 | for (int i = 0; i < MAX_NODE_SZ+1; i++) |
---|
478 | assert(thIs->child[i]->parent == thIs); |
---|
479 | for (int i = 0; i < MAX_NODE_SZ/2 + 1; i++) |
---|
480 | q->child[i] = thIs->child[i]; |
---|
481 | for (int i = MAX_NODE_SZ/2+1; i<MAX_NODE_SZ+1; i++) |
---|
482 | p->child[i-MAX_NODE_SZ/2-1] = thIs->child[i]; |
---|
483 | for (int i = 0; i < MAX_NODE_SZ/2 + 1; i++) |
---|
484 | q->child[i]->parent = q; |
---|
485 | for (int i = MAX_NODE_SZ/2+1; i < MAX_NODE_SZ+1; i++) |
---|
486 | p->child[i-MAX_NODE_SZ/2-1]->parent = p; |
---|
487 | p->update(); |
---|
488 | q->update(); |
---|
489 | |
---|
490 | if (squaredist(thIs->centre, q->centre) < squaredist(thIs->centre, p->centre)) |
---|
491 | swap(p, q); |
---|
492 | |
---|
493 | thIs->child = p->child; // now our children do not know we are their parents and believe p is their parent |
---|
494 | for (int i = 0; i < thIs->child.size(); i++) |
---|
495 | thIs->child[i]->parent = thIs; |
---|
496 | thIs->reinserted = p->reinserted; |
---|
497 | thIs->update(); |
---|
498 | delete p; |
---|
499 | |
---|
500 | if (thIs == thIs->tree->root) // root split |
---|
501 | { |
---|
502 | // since we currently are root, make new root and make us and q its children |
---|
503 | thIs->tree->newRoot(thIs->level + 1); |
---|
504 | thIs->tree->root->child.push_back(thIs); thIs->parent = thIs->tree->root; |
---|
505 | thIs->tree->root->child.push_back(q); q->parent = thIs->tree->root; |
---|
506 | thIs->tree->root->update(); |
---|
507 | } |
---|
508 | else |
---|
509 | { |
---|
510 | thIs->tree->push_front(q); |
---|
511 | } // push_front? |
---|
512 | |
---|
513 | return q; |
---|
514 | } |
---|
515 | |
---|
516 | /* Assuming we are a leaf push all leafs into our list of intersectors |
---|
517 | that are descendant of node and whoes circle intersects ours. |
---|
518 | */ |
---|
519 | void Node::search(NodePtr node) |
---|
520 | { |
---|
521 | assert(this->level == 0); |
---|
522 | int Nchild = node->child.size(); |
---|
523 | if (this->intersects(node)) { |
---|
524 | if (node->level == 0) |
---|
525 | this->intersectors.push_back(node); |
---|
526 | else |
---|
527 | for (int i=0; i<Nchild; i++) |
---|
528 | search(node->child[i]); |
---|
529 | } |
---|
530 | } |
---|
531 | |
---|
532 | /* FIXME this should not be in node.cpp and getNeighbours should not be part of the SS-tree |
---|
533 | this is mapper specific */ |
---|
534 | void findNeighbour(NodePtr thIs, NodePtr node, set<NodePtr >& neighbourList) |
---|
535 | { |
---|
536 | if (thIs->level==0) |
---|
537 | { |
---|
538 | Elt* elt1= (Elt*)(thIs->data); |
---|
539 | Elt* elt2= (Elt*)(node->data); |
---|
540 | if (isNeighbour(*elt1, *elt2)) neighbourList.insert(thIs); |
---|
541 | } |
---|
542 | else |
---|
543 | { |
---|
544 | for(int i=0; i<thIs->child.size(); i++) |
---|
545 | if (thIs->child[i]->intersects(node)) findNeighbour(thIs->child[i], node, neighbourList); |
---|
546 | } |
---|
547 | } |
---|
548 | |
---|
549 | bool Node::intersects(NodePtr node) |
---|
550 | { |
---|
551 | double d = arcdist(this->centre, node->centre); |
---|
552 | double r = this->radius; |
---|
553 | double R = node->radius; |
---|
554 | return (d < r + R + 1e-9) ? true : false; |
---|
555 | } |
---|
556 | |
---|
557 | bool Node::centreInside(Node &node) |
---|
558 | { |
---|
559 | double d = arcdist(this->centre, node.centre); |
---|
560 | double R = node.radius; |
---|
561 | return (d < R + 1e-9) ? true : false; |
---|
562 | } |
---|
563 | |
---|
564 | bool Node::isInside(Node &node) |
---|
565 | { |
---|
566 | double d = arcdist(this->centre, node.centre); |
---|
567 | double r = this->radius; |
---|
568 | double R = node.radius; |
---|
569 | return (d + r < R + 1e-9) ? true : false; |
---|
570 | } |
---|
571 | |
---|
572 | int Node::incluCheck() |
---|
573 | { |
---|
574 | if (this->level == 0) return 0; |
---|
575 | int nOutside = 0; |
---|
576 | int n = this->child.size(); // cout << "n = " << n << endl; |
---|
577 | for (int i=0; i<n; i++) |
---|
578 | { |
---|
579 | if (!this->child[i]->isInside(*this)) |
---|
580 | { |
---|
581 | cout << "Node of level " << this->level << " does not contain its " |
---|
582 | << i << "th child\n"; |
---|
583 | nOutside++; |
---|
584 | } |
---|
585 | nOutside += this->child[i]->incluCheck(); |
---|
586 | } |
---|
587 | return nOutside; |
---|
588 | } |
---|
589 | |
---|
590 | void Node::checkParent(void) |
---|
591 | { |
---|
592 | int childSize = child.size() ; |
---|
593 | |
---|
594 | for (int i = 0; i < childSize; i++) |
---|
595 | assert(child[i]->parent == this); |
---|
596 | |
---|
597 | if (level>0) for (int i = 0; i < childSize; i++) child[i]->checkParent() ; |
---|
598 | } |
---|
599 | |
---|
600 | void Node::printChildren() |
---|
601 | { |
---|
602 | cout << "level " << this->level << ", centre "; |
---|
603 | cout << "level " << this->level << ", centre " << this->centre << "\t r = " << this->radius << endl; |
---|
604 | cout << this << " p: " << this->parent << endl; |
---|
605 | int n = this->child.size(); |
---|
606 | for (int i=0; i<n; i++) |
---|
607 | { |
---|
608 | NodePtr child = this->child[i]; |
---|
609 | cout << "fils " << i << ": centre " << child->centre << "\t r = " << child->radius << endl; |
---|
610 | cout << "dist to center " << arcdist(this->centre, child->centre) << |
---|
611 | " d + R = " << arcdist(this->centre,child->centre)+child->radius << endl; |
---|
612 | } |
---|
613 | } |
---|
614 | |
---|
615 | void Node::assignRoute(std::vector<int>::iterator& rank, int level) |
---|
616 | { |
---|
617 | if (this->level==level) |
---|
618 | { |
---|
619 | route = *rank; |
---|
620 | rank++; |
---|
621 | } |
---|
622 | else |
---|
623 | { |
---|
624 | for (int i = 0; i < child.size(); i++) |
---|
625 | child[i]->assignRoute(rank, level); |
---|
626 | } |
---|
627 | } |
---|
628 | |
---|
629 | void Node::assignCircleAndPropagateUp(Coord *centres, double *radia, int level) |
---|
630 | { |
---|
631 | if (this->level == level) |
---|
632 | { |
---|
633 | // assign |
---|
634 | centre = centres[route]; |
---|
635 | radius = radia[route]; |
---|
636 | free_descendants(); // levels of sample tree beyond `level` will not be used any more |
---|
637 | child.resize(0); |
---|
638 | this->level = 0; |
---|
639 | } |
---|
640 | else |
---|
641 | { |
---|
642 | for (int i = 0; i < child.size(); i++) |
---|
643 | child[i]->assignCircleAndPropagateUp(centres, radia, level); |
---|
644 | update(); // propagate up |
---|
645 | } |
---|
646 | } |
---|
647 | |
---|
648 | /* Route node `node` within the subtree attached to us. |
---|
649 | `level` is the level one which to assign |
---|
650 | */ |
---|
651 | // Each sample node has a rank randomly assigned to it, assign `node` from full tree |
---|
652 | void Node::routeNode(NodePtr node, int level) |
---|
653 | { |
---|
654 | NodePtr closest; |
---|
655 | |
---|
656 | double distMin2 = 0; // squared distance |
---|
657 | closest = NULL; |
---|
658 | if (tree->root == this) |
---|
659 | findClosest(level, node, distMin2, closest); |
---|
660 | |
---|
661 | if (closest != NULL && tree->root == this) |
---|
662 | /* When is this point reached? |
---|
663 | if the preceeding findClosest was called and succesed to set closest |
---|
664 | findClosest sets closest if we are `level` or src is in our circle (=> belongs to child of ours) |
---|
665 | => reached if we are not `level` and node is not child of us |
---|
666 | */ |
---|
667 | node->route = closest->route; |
---|
668 | else /* find closest was not successfull or we were not root */ |
---|
669 | { |
---|
670 | if (this->level == level) |
---|
671 | node->route = this->route; |
---|
672 | else /* not yet right level => go down one more */ |
---|
673 | node->closest(this->child)->routeNode(node, level); |
---|
674 | } |
---|
675 | } |
---|
676 | |
---|
677 | void Node::routeIntersection(vector<int>& routes, NodePtr node) |
---|
678 | { |
---|
679 | if (level == 0) |
---|
680 | { |
---|
681 | routes.push_back(this->route); |
---|
682 | } |
---|
683 | else |
---|
684 | { |
---|
685 | for (int i = 0; i < child.size(); i++) { |
---|
686 | if (child[i]->intersects(node)) |
---|
687 | child[i]->routeIntersection(routes, node); |
---|
688 | } |
---|
689 | } |
---|
690 | } |
---|
691 | |
---|
692 | void Node::routingIntersecting(vector<Node> *routingList, NodePtr node) |
---|
693 | { |
---|
694 | if (level==0) |
---|
695 | { |
---|
696 | int rank = route; |
---|
697 | routingList[rank].push_back(*node); |
---|
698 | } |
---|
699 | else |
---|
700 | { |
---|
701 | for (int i = 0; i < child.size(); i++) { |
---|
702 | if (child[i]->intersects(node)) |
---|
703 | child[i]->routingIntersecting(routingList, node); |
---|
704 | } |
---|
705 | } |
---|
706 | } |
---|
707 | |
---|
708 | void Node::free_descendants() |
---|
709 | { |
---|
710 | for (int i = 0; i < child.size(); i++) |
---|
711 | { |
---|
712 | child[i]->free_descendants(); |
---|
713 | if (child[i]->level) // do not attempt to delete leafs, they are delete through leafs vector |
---|
714 | delete child[i]; |
---|
715 | } |
---|
716 | } |
---|
717 | |
---|
718 | void Node::getNodeLevel(int assignLevel, std::list<NodePtr>& NodeList) |
---|
719 | { |
---|
720 | if (level==assignLevel) NodeList.push_back(this) ; |
---|
721 | else if (level>0) for (int i = 0; i < child.size(); i++) child[i]->getNodeLevel(assignLevel,NodeList) ; |
---|
722 | return ; |
---|
723 | } |
---|
724 | |
---|
725 | bool Node::removeDeletedNodes(int assignLevel) |
---|
726 | { |
---|
727 | std::vector<NodePtr> newChild ; |
---|
728 | |
---|
729 | if (level==assignLevel+1) |
---|
730 | { |
---|
731 | bool isUpdate=false ; |
---|
732 | for (int i = 0; i < child.size(); i++) |
---|
733 | { |
---|
734 | if (child[i]->toDelete) |
---|
735 | { |
---|
736 | isUpdate=true ; |
---|
737 | for (int j = 0; j < child[i]->child.size(); j++) tree->push_back(child[i]->child[j]) ; |
---|
738 | tree->decreaseLevelSize(assignLevel) ; |
---|
739 | delete child[i] ; |
---|
740 | } |
---|
741 | else newChild.push_back(child[i]) ; |
---|
742 | } |
---|
743 | |
---|
744 | if (isUpdate) |
---|
745 | { |
---|
746 | child=newChild ; |
---|
747 | update() ; |
---|
748 | return true ; |
---|
749 | } |
---|
750 | else return false ; |
---|
751 | } |
---|
752 | else |
---|
753 | { |
---|
754 | bool isUpdate=false ; |
---|
755 | for (int i = 0; i < child.size(); i++) isUpdate |= child[i]->removeDeletedNodes(assignLevel) ; |
---|
756 | if (isUpdate) update() ; |
---|
757 | return isUpdate ; |
---|
758 | } |
---|
759 | } |
---|
760 | |
---|
761 | } |
---|