All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
QuadNodeCartesianEuclid.h
Go to the documentation of this file.
1 /*
2  * QuadNodePolarEuclid.h
3  *
4  * Created on: 21.05.2014
5  * Author: Moritz v. Looz (moritz.looz-corswarem@kit.edu)
6  *
7  * Note: This is similar enough to QuadNode.h that one could merge these two classes.
8  */
9 
10 #ifndef QUADNODECARTESIANEUCLID_H_
11 #define QUADNODECARTESIANEUCLID_H_
12 
13 #include <vector>
14 #include <algorithm>
15 #include <functional>
16 #include <assert.h>
17 #include "../../auxiliary/Log.h"
18 #include "../../geometric/HyperbolicSpace.h"
19 
20 using std::vector;
21 using std::min;
22 using std::max;
23 using std::cos;
24 
25 namespace NetworKit {
26 
27 template <class T>
29  friend class QuadTreeGTest;
30 private:
31  Point<double> minPoint;
32  Point<double> maxPoint;
33  count dimension;
34  unsigned capacity;
35  static const unsigned coarsenLimit = 4;
36  count subTreeSize;
37  std::vector<T> content;
38  std::vector<Point<double> > positions;
39  bool isLeaf;
40  bool splitTheoretical;
41  index ID;
42  double lowerBoundR;
43 
44 public:
45  std::vector<QuadNodeCartesianEuclid> children;
46 
57  QuadNodeCartesianEuclid(Point<double> lower = Point<double>({0.0, 0.0}), Point<double> upper = Point<double>({1.0, 1.0}), unsigned capacity = 1000, bool splitTheoretical = false) {
58  this->minPoint = lower;
59  this->maxPoint = upper;
60  this->dimension = minPoint.getDimensions();
61  assert(maxPoint.getDimensions() == dimension);
62  this->capacity = capacity;
63  this->splitTheoretical = splitTheoretical;
64  this->ID = 0;
65  isLeaf = true;
66  subTreeSize = 0;
67  }
68 
69  void split() {
70  assert(isLeaf);
71  assert(children.size() == 0);
72  vector<double> middle(dimension);
73  if (splitTheoretical) {
74  //Euclidean space is distributed equally
75  for (index d = 0; d < dimension; d++) {
76  middle[d] = (minPoint[d] + maxPoint[d]) / 2;
77  }
78  } else {
79  //median of points
80  const count numPoints = positions.size();
81  assert(numPoints > 0);//otherwise, why split?
82  vector<vector<double> > sorted(dimension);
83  for (index d = 0; d < dimension; d++) {
84  sorted[d].resize(numPoints);
85  for (index i = 0; i < numPoints; i++) {
86  sorted[d][i] = positions[i][d];
87  }
88  std::sort(sorted[d].begin(), sorted[d].end());
89  middle[d] = sorted[d][numPoints/2];//this will crash if no points are there!
90  assert(middle[d] <= maxPoint[d]);
91  assert(middle[d] >= minPoint[d]);
92  }
93  }
94  count childCount = pow(2,dimension);
95  for (index i = 0; i < childCount; i++) {
96  vector<double> lowerValues(dimension);
97  vector<double> upperValues(dimension);
98  index bitCopy = i;
99  for (index d = 0; d < dimension; d++) {
100  if (bitCopy & 1) {
101  lowerValues[d] = middle[d];
102  upperValues[d] = maxPoint[d];
103  } else {
104  lowerValues[d] = minPoint[d];
105  upperValues[d] = middle[d];
106  }
107  bitCopy = bitCopy >> 1;
108  }
109  QuadNodeCartesianEuclid child(Point<double>(lowerValues), Point<double>(upperValues), capacity, splitTheoretical);
110  assert(child.isLeaf);
111  children.push_back(child);
112  }
113  isLeaf = false;
114  }
115 
123  void addContent(T input, Point<double> pos) {
124  assert(content.size() == positions.size());
125  assert(this->responsible(pos));
126  if (isLeaf) {
127  if (content.size() + 1 < capacity) {
128  content.push_back(input);
129  positions.push_back(pos);
130  } else {
131  split();
132 
133  for (index i = 0; i < content.size(); i++) {
134  this->addContent(content[i], positions[i]);
135  }
136  assert(subTreeSize == content.size());//we have added everything twice
137  subTreeSize = content.size();
138  content.clear();
139  positions.clear();
140  this->addContent(input, pos);
141  }
142  }
143  else {
144  assert(children.size() > 0);
145  bool foundResponsibleChild = false;
146  for (index i = 0; i < children.size(); i++) {
147  if (children[i].responsible(pos)) {
148  foundResponsibleChild = true;
149  children[i].addContent(input, pos);
150  break;
151  }
152  }
153  assert(foundResponsibleChild);
154  subTreeSize++;
155  }
156  }
157 
166  bool removeContent(T input, Point<double> pos) {
167  if (!responsible(pos)) return false;
168  if (isLeaf) {
169  index i = 0;
170  for (; i < content.size(); i++) {
171  if (content[i] == input) break;
172  }
173  if (i < content.size()) {
174  assert(positions[i].distance(pos) == 0);
175  //remove element
176  content.erase(content.begin()+i);
177  positions.erase(positions.begin()+i);
178  return true;
179  } else {
180  return false;
181  }
182  }
183  else {
184  bool removed = false;
185  bool allLeaves = true;
186  assert(children.size() > 0);
187  for (index i = 0; i < children.size(); i++) {
188  if (!children[i].isLeaf) allLeaves = false;
189  if (children[i].removeContent(input, pos)) {
190  assert(!removed);
191  removed = true;
192  }
193  }
194  if (removed) subTreeSize--;
195  //coarsen?
196  if (removed && allLeaves && size() < coarsenLimit) {
197  //coarsen!!
198  //why not assert empty containers and then insert directly?
199  vector<T> allContent;
200  vector<Point<double> > allPositions;
201  for (index i = 0; i < children.size(); i++) {
202  allContent.insert(allContent.end(), children[i].content.begin(), children[i].content.end());
203  allPositions.insert(allPositions.end(), children[i].positions.begin(), children[i].positions.end());
204  }
205  assert(allContent.size() == allPositions.size());
206  children.clear();
207  content.swap(allContent);
208  positions.swap(allPositions);
209  isLeaf = true;
210  }
211 
212  return removed;
213  }
214  }
215 
216 
225  bool outOfReach(Point<double> query, double radius) const {
226  return EuclideanDistances(query).first > radius;
227  }
228 
232  std::pair<double, double> EuclideanDistances(Point<double> query) const {
237  double maxDistance = 0;
238  double minDistance = std::numeric_limits<double>::max();
239  //Point<double> minCopy(minPoint);
240  //Point<double> maxCopy(minPoint);
241 
242  if (responsible(query)) minDistance = 0;
243 
244  auto updateMinMax = [&minDistance, &maxDistance, query](Point<double> pos){
245  double extremalValue = pos.distance(query);
246  maxDistance = std::max(extremalValue, maxDistance);
247  minDistance = std::min(minDistance, extremalValue);
248  };
249 
250  vector<double> closestValues(dimension);
251  vector<double> farthestValues(dimension);
252 
253  for (index d = 0; d < dimension; d++) {
254  if (std::abs(query[d] - minPoint.at(d)) < std::abs(query[d] - maxPoint.at(d))) {
255  closestValues[d] = minPoint.at(d);
256  farthestValues[d] = maxPoint.at(d);
257  } else {
258  farthestValues[d] = minPoint.at(d);
259  closestValues[d] = maxPoint.at(d);
260  }
261  if (query[d] >= minPoint.at(d) && query[d] <= maxPoint.at(d)) {
262  closestValues[d] = query[d];
263  }
264  }
265  updateMinMax(Point<double>(closestValues));
266  updateMinMax(Point<double>(farthestValues));
267 
268  assert(minDistance < query.length() + maxPoint.length());
269  assert(minDistance < maxDistance);
270  return std::pair<double, double>(minDistance, maxDistance);
271  }
272 
273 
282  bool responsible(Point<double> pos) const {
283  for (index d = 0; d < dimension; d++) {
284  if (pos[d] < minPoint.at(d) || pos[d] >= maxPoint.at(d)) return false;
285  }
286  return true;
287  }
288 
294  std::vector<T> getElements() const {
295  if (isLeaf) {
296  return content;
297  } else {
298  assert(content.size() == 0);
299  assert(positions.size() == 0);
300  vector<T> result;
301  for (index i = 0; i < children.size(); i++) {
302  std::vector<T> subresult = children[i].getElements();
303  result.insert(result.end(), subresult.begin(), subresult.end());
304  }
305  return result;
306  }
307  }
308 
309  void getCoordinates(vector<Point<double> > &pointContainer) const {
310  if (isLeaf) {
311  pointContainer.insert(pointContainer.end(), positions.begin(), positions.end());
312  }
313  else {
314  assert(content.size() == 0);
315  assert(positions.size() == 0);
316  for (index i = 0; i < children.size(); i++) {
317  children[i].getCoordinates(pointContainer);
318  }
319  }
320  }
321 
339  void getElementsInEuclideanCircle(Point<double> center, double radius, vector<T> &result) const {
340  if (outOfReach(center, radius)) {
341  return;
342  }
343 
344  if (isLeaf) {
345  const double rsq = radius*radius;
346  const count cSize = content.size();
347 
348  for (index i=0; i < cSize; i++) {
349  if (positions[i].squaredDistance(center) < rsq) {
350  result.push_back(content[i]);
351  }
352  }
353  } else {
354  for (index i = 0; i < children.size(); i++) {
355  children[i].getElementsInEuclideanCircle(center, radius, result);
356  }
357  }
358  }
359 
360  count getElementsProbabilistically(Point<double> euQuery, std::function<double(double)> prob, vector<T> &result) const {
361  TRACE("Getting Euclidean distances");
362  auto distancePair = EuclideanDistances(euQuery);
363  double probUB = prob(distancePair.first);
364  double probLB = prob(distancePair.second);
365  assert(probLB <= probUB);
366  if (probUB > 0.5) probUB = 1;
367  if (probUB == 0) return 0;
368  //TODO: return whole if probLB == 1
369  double probdenom = std::log(1-probUB);
370  if (probdenom == 0) return 0;//there is a very small probability, but we cannot process it.
371  TRACE("probUB: ", probUB, ", probdenom: ", probdenom);
372 
373  count expectedNeighbours = probUB*size();
374  count candidatesTested = 0;
375  count incomingNeighbours = result.size();
376  count ownsize = size();
377 
378 
379  if (isLeaf) {
380  const count lsize = content.size();
381  TRACE("Leaf of size ", lsize);
382  for (index i = 0; i < lsize; i++) {
383  //jump!
384  if (probUB < 1) {
385  double random = Aux::Random::real();
386  double delta = std::log(random) / probdenom;
387  assert(delta >= 0);
388  i += delta;
389  if (i >= lsize) break;
390  TRACE("Jumped with delta ", delta, " arrived at ", i);
391  }
392  assert(i >= 0);
393 
394  //see where we've arrived
395  candidatesTested++;
396  double distance = positions[i].distance(euQuery);
397  assert(distance >= distancePair.first);//TODO: These should not fail!
398  assert(distance <= distancePair.second);
399  double q = prob(distance);
400  q = q / probUB; //since the candidate was selected by the jumping process, we have to adjust the probabilities
401  assert(q <= 1);
402 
403  //accept?
404  double acc = Aux::Random::real();
405  if (acc < q) {
406  TRACE("Accepted node ", i, " with probability ", q, ".");
407  result.push_back(content[i]);
408  }
409  }
410  } else {
411  if (expectedNeighbours < 4 || probUB < 1/1000) {//select candidates directly instead of calling recursively
412  TRACE("probUB = ", probUB, ", switching to direct candidate selection.");
413  assert(probUB < 1);
414  const count stsize = size();
415  for (index i = 0; i < stsize; i++) {
416  double delta = std::log(Aux::Random::real()) / probdenom;
417  assert(delta >= 0);
418  i += delta;
419  TRACE("Jumped with delta ", delta, " arrived at ", i, ". Calling maybeGetKthElement.");
420  if (i < size()) maybeGetKthElement(probUB, euQuery, prob, i, result);//this could be optimized. As of now, the offset is subtracted separately for each point
421  else break;
422  candidatesTested++;
423  }
424  } else {//carry on as normal
425  for (index i = 0; i < children.size(); i++) {
426  TRACE("Recursively calling child ", i);
427  candidatesTested += children[i].getElementsProbabilistically(euQuery, prob, result);
428  }
429  }
430  }
431  count finalNeighbours = result.size();
432  if (probLB == 1) assert(finalNeighbours == incomingNeighbours + ownsize);
433  return candidatesTested;
434  }
435 
436 
437  void maybeGetKthElement(double upperBound, Point<double> euQuery, std::function<double(double)> prob, index k, vector<T> &circleDenizens) const {
438  TRACE("Maybe get element ", k, " with upper Bound ", upperBound);
439  assert(k < size());
440  if (isLeaf) {
441  double acceptance = prob(euQuery.distance(positions[k]))/upperBound;
442  TRACE("Is leaf, accept with ", acceptance);
443  if (Aux::Random::real() < acceptance) circleDenizens.push_back(content[k]);
444  } else {
445  TRACE("Call recursively.");
446  index offset = 0;
447  for (index i = 0; i < children.size(); i++) {
448  count childsize = children[i].size();
449  if (k - offset < childsize) {
450  children[i].maybeGetKthElement(upperBound, euQuery, prob, k - offset, circleDenizens);
451  break;
452  }
453  offset += childsize;
454  }
455  }
456  }
457 
462  void trim() {
463  content.shrink_to_fit();
464  positions.shrink_to_fit();
465  if (!isLeaf) {
466  for (index i = 0; i < children.size(); i++) {
467  children[i].trim();
468  }
469  }
470  }
471 
475  count size() const {
476  return isLeaf ? content.size() : subTreeSize;
477  }
478 
479  void recount() {
480  subTreeSize = 0;
481  for (index i = 0; i < children.size(); i++) {
482  children[i].recount();
483  subTreeSize += children[i].size();
484  }
485  }
486 
490  count height() const {
491  count result = 1;//if leaf node, the children loop will not execute
492  for (auto child : children) result = std::max(result, child.height()+1);
493  return result;
494  }
495 
499  count countLeaves() const {
500  if (isLeaf) return 1;
501  count result = 0;
502  for (index i = 0; i < children.size(); i++) {
503  result += children[i].countLeaves();
504  }
505  return result;
506  }
507 
508  index getID() const {
509  return ID;
510  }
511 
513  index result = nextID;
514  assert(children.size() == pow(2,dimension) || children.size() == 0);
515  for (int i = 0; i < children.size(); i++) {
516  result = children[i].indexSubtree(result);
517  }
518  this->ID = result;
519  return result+1;
520  }
521 
523  if (!responsible(pos)) return NetworKit::none;
524  if (isLeaf) return getID();
525  else {
526  for (int i = 0; i < children.size(); i++) {
527  index childresult = children[i].getCellID(pos);
528  if (childresult != NetworKit::none) return childresult;
529  }
530  throw std::runtime_error("No responsible child node found even though this node is responsible.");
531  }
532  }
533 
535  if (isLeaf) return getID();
536  else {
537  index result = -1;
538  for (int i = 0; i < children.size(); i++) {
539  result = std::max(children[i].getMaxIDInSubtree(), result);
540  }
541  return std::max(result, getID());
542  }
543  }
544 
545  count reindex(count offset) {
546  if (isLeaf)
547  {
548  #pragma omp task
549  {
550  index p = offset;
551  std::generate(content.begin(), content.end(), [&p](){return p++;});
552  }
553  offset += size();
554  } else {
555  for (int i = 0; i < children.size(); i++) {
556  offset = children[i].reindex(offset);
557  }
558  }
559  return offset;
560  }
561 };
562 }
563 
564 #endif /* QUADNODE_H_ */
void maybeGetKthElement(double upperBound, Point< double > euQuery, std::function< double(double)> prob, index k, vector< T > &circleDenizens) const
Definition: QuadNodeCartesianEuclid.h:437
index getCellID(Point< double > pos) const
Definition: QuadNodeCartesianEuclid.h:522
QuadNodeCartesianEuclid(Point< double > lower=Point< double >({0.0, 0.0}), Point< double > upper=Point< double >({1.0, 1.0}), unsigned capacity=1000, bool splitTheoretical=false)
Construct a QuadNode for polar coordinates.
Definition: QuadNodeCartesianEuclid.h:57
double real()
Definition: Random.cpp:77
index getMaxIDInSubtree() const
Definition: QuadNodeCartesianEuclid.h:534
std::pair< double, double > EuclideanDistances(Point< double > query) const
Definition: QuadNodeCartesianEuclid.h:232
void getElementsInEuclideanCircle(Point< double > center, double radius, vector< T > &result) const
Main query method, get points lying in a Euclidean circle around the center point.
Definition: QuadNodeCartesianEuclid.h:339
void log(const Location &loc, LogLevel p, const std::string msg)
Definition: Log.cpp:202
T length() const
Definition: Point.h:99
bool responsible(Point< double > pos) const
Does the point at (angle, r) fall inside the region managed by this QuadNode?
Definition: QuadNodeCartesianEuclid.h:282
uint64_t index
Typedefs.
Definition: Globals.h:20
count reindex(count offset)
Definition: QuadNodeCartesianEuclid.h:545
void trim()
Shrink all vectors in this subtree to fit the content.
Definition: QuadNodeCartesianEuclid.h:462
friend class QuadTreeGTest
Definition: QuadNodeCartesianEuclid.h:29
Definition: QuadNodeCartesianEuclid.h:28
std::vector< T > getElements() const
Get all Elements in this QuadNode or a descendant of it.
Definition: QuadNodeCartesianEuclid.h:294
std::vector< QuadNodeCartesianEuclid > children
Definition: QuadNodeCartesianEuclid.h:45
void recount()
Definition: QuadNodeCartesianEuclid.h:479
T & at(const index i)
Definition: Point.h:217
T distance(const Point< T > &p) const
Definition: Point.h:130
index getID() const
Definition: QuadNodeCartesianEuclid.h:508
uint64_t count
Definition: Globals.h:21
void getCoordinates(vector< Point< double > > &pointContainer) const
Definition: QuadNodeCartesianEuclid.h:309
constexpr index none
Constants.
Definition: Globals.h:28
bool removeContent(T input, Point< double > pos)
Remove content at coordinate pos.
Definition: QuadNodeCartesianEuclid.h:166
count getElementsProbabilistically(Point< double > euQuery, std::function< double(double)> prob, vector< T > &result) const
Definition: QuadNodeCartesianEuclid.h:360
count countLeaves() const
Leaf cells in the subtree hanging from this QuadNode.
Definition: QuadNodeCartesianEuclid.h:499
float distance
Definition: PubWebGenerator.h:25
count height() const
Height of subtree hanging from this QuadNode.
Definition: QuadNodeCartesianEuclid.h:490
void split()
Definition: QuadNodeCartesianEuclid.h:69
index indexSubtree(index nextID)
Definition: QuadNodeCartesianEuclid.h:512
bool outOfReach(Point< double > query, double radius) const
Check whether the region managed by this node lies outside of an Euclidean circle.
Definition: QuadNodeCartesianEuclid.h:225
#define TRACE(...)
Definition: Log.h:92
void addContent(T input, Point< double > pos)
Add a point at polar coordinates (angle, R) with content input.
Definition: QuadNodeCartesianEuclid.h:123
count getDimensions() const
Definition: Point.h:50
count size() const
Number of points lying in the region managed by this QuadNode.
Definition: QuadNodeCartesianEuclid.h:475