Tapkee
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
quadtree.hpp
Go to the documentation of this file.
1 
32 #include <math.h>
33 #include <float.h>
34 #include <stdlib.h>
35 #include <stdio.h>
36 #include <algorithm>
37 
38 #ifndef QUADTREE_H
39 #define QUADTREE_H
40 
41 namespace tsne
42 {
43 
44 class Cell {
45 
46 public:
47 
48  double x;
49  double y;
50  double hw;
51  double hh;
52 
53  bool containsPoint(double point[])
54  {
55  if(x - hw > point[0]) return false;
56  if(x + hw < point[0]) return false;
57  if(y - hh > point[1]) return false;
58  if(y + hh < point[1]) return false;
59  return true;
60  }
61 
62 };
63 
64 
65 class QuadTree
66 {
67 
68  // Fixed constants
69  static const int QT_NO_DIMS = 2;
70  static const int QT_NODE_CAPACITY = 1;
71 
72  // A buffer we use when doing force computations
73  double buff[QT_NO_DIMS];
74 
75  // Properties of this node in the tree
77  bool is_leaf;
78  int size;
79  int cum_size;
80 
81  // Axis-aligned bounding box stored as a center with half-dimensions to represent the boundaries of this quad tree
83 
84  // Indices in this quad tree node, corresponding center-of-mass, and list of all children
85  double* data;
88 
89  // Children
94 
95 public:
96 
97  // Default constructor for quadtree -- build tree, too!
98  QuadTree(double* inp_data, int N) :
99  parent(NULL), is_leaf(false), size(0), cum_size(0), boundary(), data(NULL),
100  northWest(NULL), northEast(NULL), southWest(NULL), southEast(NULL)
101  {
102  // Compute mean, width, and height of current map (boundaries of quadtree)
103  double* mean_Y = new double[QT_NO_DIMS]; for(int d = 0; d < QT_NO_DIMS; d++) mean_Y[d] = .0;
104  double* min_Y = new double[QT_NO_DIMS]; for(int d = 0; d < QT_NO_DIMS; d++) min_Y[d] = DBL_MAX;
105  double* max_Y = new double[QT_NO_DIMS]; for(int d = 0; d < QT_NO_DIMS; d++) max_Y[d] = -DBL_MAX;
106  for(int n = 0; n < N; n++) {
107  for(int d = 0; d < QT_NO_DIMS; d++) {
108  mean_Y[d] += inp_data[n * QT_NO_DIMS + d];
109  if(inp_data[n * QT_NO_DIMS + d] < min_Y[d]) min_Y[d] = inp_data[n * QT_NO_DIMS + d];
110  if(inp_data[n * QT_NO_DIMS + d] > max_Y[d]) max_Y[d] = inp_data[n * QT_NO_DIMS + d];
111  }
112  }
113  for(int d = 0; d < QT_NO_DIMS; d++) mean_Y[d] /= (double) N;
114 
115  // Construct quadtree
116  init(NULL, inp_data, mean_Y[0], mean_Y[1], std::max(max_Y[0] - mean_Y[0], mean_Y[0] - min_Y[0]) + 1e-5,
117  std::max(max_Y[1] - mean_Y[1], mean_Y[1] - min_Y[1]) + 1e-5);
118  fill(N);
119  delete[] mean_Y; delete[] max_Y; delete[] min_Y;
120  }
121 
122  // Constructor for quadtree with particular size and parent -- build the tree, too!
123  QuadTree(double* inp_data, double inp_x, double inp_y, double inp_hw, double inp_hh) :
124  parent(NULL), is_leaf(false), size(0), cum_size(0), boundary(), data(NULL),
125  northWest(NULL), northEast(NULL), southWest(NULL), southEast(NULL)
126  {
127  init(NULL, inp_data, inp_x, inp_y, inp_hw, inp_hh);
128  }
129 
130  // Constructor for quadtree with particular size and parent -- build the tree, too!
131  QuadTree(double* inp_data, int N, double inp_x, double inp_y, double inp_hw, double inp_hh) :
132  parent(NULL), is_leaf(false), size(0), cum_size(0), boundary(), data(NULL),
133  northWest(NULL), northEast(NULL), southWest(NULL), southEast(NULL)
134  {
135  init(NULL, inp_data, inp_x, inp_y, inp_hw, inp_hh);
136  fill(N);
137  }
138 
139  // Constructor for quadtree with particular size (do not fill the tree)
140  QuadTree(QuadTree* inp_parent, double* inp_data, int N, double inp_x, double inp_y, double inp_hw, double inp_hh) :
141  parent(NULL), is_leaf(false), size(0), cum_size(0), boundary(), data(NULL),
142  northWest(NULL), northEast(NULL), southWest(NULL), southEast(NULL)
143  {
144  init(inp_parent, inp_data, inp_x, inp_y, inp_hw, inp_hh);
145  fill(N);
146  }
147 
148  // Constructor for quadtree with particular size and parent (do not fill the tree)
149  QuadTree(QuadTree* inp_parent, double* inp_data, double inp_x, double inp_y, double inp_hw, double inp_hh) :
150  parent(NULL), is_leaf(false), size(0), cum_size(0), boundary(), data(NULL),
151  northWest(NULL), northEast(NULL), southWest(NULL), southEast(NULL)
152  {
153  init(inp_parent, inp_data, inp_x, inp_y, inp_hw, inp_hh);
154  }
155 
156  // Destructor for quadtree
158  {
159  delete northWest;
160  delete northEast;
161  delete southWest;
162  delete southEast;
163  }
164 
165  void setData(double* inp_data)
166  {
167  data = inp_data;
168  }
169 
171  {
172  return parent;
173  }
174 
175  //void construct(Cell boundary);
176 
177  // Insert a point into the QuadTree
178  bool insert(int new_index)
179  {
180  // Ignore objects which do not belong in this quad tree
181  double* point = data + new_index * QT_NO_DIMS;
182  if(!boundary.containsPoint(point))
183  return false;
184 
185  // Online update of cumulative size and center-of-mass
186  cum_size++;
187  double mult1 = (double) (cum_size - 1) / (double) cum_size;
188  double mult2 = 1.0 / (double) cum_size;
189  for(int d = 0; d < QT_NO_DIMS; d++) center_of_mass[d] *= mult1;
190  for(int d = 0; d < QT_NO_DIMS; d++) center_of_mass[d] += mult2 * point[d];
191 
192  // If there is space in this quad tree and it is a leaf, add the object here
193  if(is_leaf && size < QT_NODE_CAPACITY) {
194  index[size] = new_index;
195  size++;
196  return true;
197  }
198 
199  // Don't add duplicates for now (this is not very nice)
200  bool any_duplicate = false;
201  for(int n = 0; n < size; n++) {
202  bool duplicate = true;
203  for(int d = 0; d < QT_NO_DIMS; d++) {
204  if(point[d] != data[index[n] * QT_NO_DIMS + d]) { duplicate = false; break; }
205  }
206  any_duplicate = any_duplicate | duplicate;
207  }
208  if(any_duplicate) return true;
209 
210  // Otherwise, we need to subdivide the current cell
211  if(is_leaf) subdivide();
212 
213  // Find out where the point can be inserted
214  if(northWest->insert(new_index)) return true;
215  if(northEast->insert(new_index)) return true;
216  if(southWest->insert(new_index)) return true;
217  if(southEast->insert(new_index)) return true;
218 
219  // Otherwise, the point cannot be inserted (this should never happen)
220  return false;
221  }
222 
223  // Create four children which fully divide this cell into four quads of equal area
224  void subdivide()
225  {
226  // Create four children
227  northWest = new QuadTree(this, data, boundary.x - .5 * boundary.hw, boundary.y - .5 * boundary.hh, .5 * boundary.hw, .5 * boundary.hh);
228  northEast = new QuadTree(this, data, boundary.x + .5 * boundary.hw, boundary.y - .5 * boundary.hh, .5 * boundary.hw, .5 * boundary.hh);
229  southWest = new QuadTree(this, data, boundary.x - .5 * boundary.hw, boundary.y + .5 * boundary.hh, .5 * boundary.hw, .5 * boundary.hh);
230  southEast = new QuadTree(this, data, boundary.x + .5 * boundary.hw, boundary.y + .5 * boundary.hh, .5 * boundary.hw, .5 * boundary.hh);
231 
232  // Move existing points to correct children
233  for(int i = 0; i < size; i++) {
234  bool success = false;
235  if(!success) success = northWest->insert(index[i]);
236  if(!success) success = northEast->insert(index[i]);
237  if(!success) success = southWest->insert(index[i]);
238  if(!success) success = southEast->insert(index[i]);
239  index[i] = -1;
240  }
241 
242  // Empty parent node
243  size = 0;
244  is_leaf = false;
245  }
246 
247  // Checks whether the specified tree is correct
248  bool isCorrect()
249  {
250  for(int n = 0; n < size; n++) {
251  double* point = data + index[n] * QT_NO_DIMS;
252  if(!boundary.containsPoint(point)) return false;
253  }
254  if(!is_leaf) return northWest->isCorrect() &&
255  northEast->isCorrect() &&
256  southWest->isCorrect() &&
257  southEast->isCorrect();
258  else return true;
259  }
260 
261  // Rebuilds a possibly incorrect tree (LAURENS: This function is not tested yet!)
262  void rebuildTree()
263  {
264  for(int n = 0; n < size; n++) {
265  // Check whether point is erroneous
266  double* point = data + index[n] * QT_NO_DIMS;
267  if(!boundary.containsPoint(point)) {
268 
269  // Remove erroneous point
270  int rem_index = index[n];
271  for(int m = n + 1; m < size; m++) index[m - 1] = index[m];
272  index[size - 1] = -1;
273  size--;
274 
275  // Update center-of-mass and counter in all parents
276  bool done = false;
277  QuadTree* node = this;
278  while(!done) {
279  for(int d = 0; d < QT_NO_DIMS; d++) {
280  node->center_of_mass[d] = ((double) node->cum_size * node->center_of_mass[d] - point[d]) / (double) (node->cum_size - 1);
281  }
282  node->cum_size--;
283  if(node->getParent() == NULL) done = true;
284  else node = node->getParent();
285  }
286 
287  // Reinsert point in the root tree
288  node->insert(rem_index);
289  }
290  }
291 
292  // Rebuild lower parts of the tree
297  }
298 
299  // Build a list of all indices in quadtree
300  void getAllIndices(int* indices)
301  {
302  getAllIndices(indices, 0);
303  }
304 
305  int getDepth()
306  {
307  if(is_leaf) return 1;
308  return 1 + std::max(std::max(northWest->getDepth(),
309  northEast->getDepth()),
310  std::max(southWest->getDepth(),
311  southEast->getDepth()));
312  }
313 
314  // Compute non-edge forces using Barnes-Hut algorithm
315  void computeNonEdgeForces(int point_index, double theta, double neg_f[], double* sum_Q)
316  {
317 
318  // Make sure that we spend no time on empty nodes or self-interactions
319  if(cum_size == 0 || (is_leaf && size == 1 && index[0] == point_index)) return;
320 
321  // Compute distance between point and center-of-mass
322  double D = .0;
323  int ind = point_index * QT_NO_DIMS;
324  for(int d = 0; d < QT_NO_DIMS; d++) buff[d] = data[ind + d];
325  for(int d = 0; d < QT_NO_DIMS; d++) buff[d] -= center_of_mass[d];
326  for(int d = 0; d < QT_NO_DIMS; d++) D += buff[d] * buff[d];
327 
328  // Check whether we can use this node as a "summary"
329  if(is_leaf || std::max(boundary.hh, boundary.hw)/sqrt(D) < theta) {
330 
331  // Compute and add t-SNE force between point and current node
332  double Q = 1.0 / (1.0 + D);
333  *sum_Q += cum_size * Q;
334  double mult = cum_size * Q * Q;
335  for(int d = 0; d < QT_NO_DIMS; d++) neg_f[d] += mult * buff[d];
336  }
337  else {
338 
339  // Recursively apply Barnes-Hut to children
340  northWest->computeNonEdgeForces(point_index, theta, neg_f, sum_Q);
341  northEast->computeNonEdgeForces(point_index, theta, neg_f, sum_Q);
342  southWest->computeNonEdgeForces(point_index, theta, neg_f, sum_Q);
343  southEast->computeNonEdgeForces(point_index, theta, neg_f, sum_Q);
344  }
345  }
346 
347  // Computes edge forces
348  void computeEdgeForces(int* row_P, int* col_P, double* val_P, int N, double* pos_f)
349  {
350  // Loop over all edges in the graph
351  int ind1, ind2;
352  double D;
353  for(int n = 0; n < N; n++) {
354  ind1 = n * QT_NO_DIMS;
355  for(int i = row_P[n]; i < row_P[n + 1]; i++) {
356 
357  // Compute pairwise distance and Q-value
358  D = .0;
359  ind2 = col_P[i] * QT_NO_DIMS;
360  for(int d = 0; d < QT_NO_DIMS; d++) buff[d] = data[ind1 + d];
361  for(int d = 0; d < QT_NO_DIMS; d++) buff[d] -= data[ind2 + d];
362  for(int d = 0; d < QT_NO_DIMS; d++) D += buff[d] * buff[d];
363  D = val_P[i] / (1.0 + D);
364 
365  // Sum positive force
366  for(int d = 0; d < QT_NO_DIMS; d++) pos_f[ind1 + d] += D * buff[d];
367  }
368  }
369  }
370 
371  // Print out tree
372  void print()
373  {
374  if(cum_size == 0) {
375  printf("Empty node\n");
376  return;
377  }
378 
379  if(is_leaf) {
380  printf("Leaf node; data = [");
381  for(int i = 0; i < size; i++) {
382  double* point = data + index[i] * QT_NO_DIMS;
383  for(int d = 0; d < QT_NO_DIMS; d++) printf("%f, ", point[d]);
384  printf(" (index = %d)", index[i]);
385  if(i < size - 1) printf("\n");
386  else printf("]\n");
387  }
388  }
389  else {
390  printf("Intersection node with center-of-mass = [");
391  for(int d = 0; d < QT_NO_DIMS; d++) printf("%f, ", center_of_mass[d]);
392  printf("]; children are:\n");
393  northEast->print();
394  northWest->print();
395  southEast->print();
396  southWest->print();
397  }
398  }
399 
400 private:
401 
402  QuadTree(const QuadTree&);
403  QuadTree& operator=(const QuadTree&);
404 
405  void init(QuadTree* inp_parent, double* inp_data, double inp_x, double inp_y, double inp_hw, double inp_hh)
406  {
407  parent = inp_parent;
408  data = inp_data;
409  is_leaf = true;
410  size = 0;
411  cum_size = 0;
412  boundary.x = inp_x;
413  boundary.y = inp_y;
414  boundary.hw = inp_hw;
415  boundary.hh = inp_hh;
416  northWest = NULL;
417  northEast = NULL;
418  southWest = NULL;
419  southEast = NULL;
420  for(int i = 0; i < QT_NO_DIMS; i++) center_of_mass[i] = .0;
421  }
422 
423  // Build quadtree on dataset
424  void fill(int N)
425  {
426  for(int i = 0; i < N; i++) insert(i);
427  }
428 
429  // Build a list of all indices in quadtree
430  int getAllIndices(int* indices, int loc)
431  {
432 
433  // Gather indices in current quadrant
434  for(int i = 0; i < size; i++) indices[loc + i] = index[i];
435  loc += size;
436 
437  // Gather indices in children
438  if(!is_leaf) {
439  loc = northWest->getAllIndices(indices, loc);
440  loc = northEast->getAllIndices(indices, loc);
441  loc = southWest->getAllIndices(indices, loc);
442  loc = southEast->getAllIndices(indices, loc);
443  }
444  return loc;
445  }
446 
447  //bool isChild(int test_index, int start, int end);
448 };
449 
450 }
451 
452 #endif
QuadTree * southWest
Definition: quadtree.hpp:92
QuadTree * parent
Definition: quadtree.hpp:76
double y
Definition: quadtree.hpp:49
QuadTree * northWest
Definition: quadtree.hpp:90
double hw
Definition: quadtree.hpp:50
void rebuildTree()
Definition: quadtree.hpp:262
void getAllIndices(int *indices)
Definition: quadtree.hpp:300
QuadTree * northEast
Definition: quadtree.hpp:91
bool containsPoint(double point[])
Definition: quadtree.hpp:53
double center_of_mass[QT_NO_DIMS]
Definition: quadtree.hpp:86
void setData(double *inp_data)
Definition: quadtree.hpp:165
bool isCorrect()
Definition: quadtree.hpp:248
static const int QT_NODE_CAPACITY
Definition: quadtree.hpp:70
double hh
Definition: quadtree.hpp:51
void computeEdgeForces(int *row_P, int *col_P, double *val_P, int N, double *pos_f)
Definition: quadtree.hpp:348
QuadTree(double *inp_data, int N, double inp_x, double inp_y, double inp_hw, double inp_hh)
Definition: quadtree.hpp:131
QuadTree * getParent()
Definition: quadtree.hpp:170
QuadTree(QuadTree *inp_parent, double *inp_data, double inp_x, double inp_y, double inp_hw, double inp_hh)
Definition: quadtree.hpp:149
QuadTree * southEast
Definition: quadtree.hpp:93
int index[QT_NODE_CAPACITY]
Definition: quadtree.hpp:87
static const int QT_NO_DIMS
Definition: quadtree.hpp:69
bool insert(int new_index)
Definition: quadtree.hpp:178
int getAllIndices(int *indices, int loc)
Definition: quadtree.hpp:430
QuadTree(QuadTree *inp_parent, double *inp_data, int N, double inp_x, double inp_y, double inp_hw, double inp_hh)
Definition: quadtree.hpp:140
double buff[QT_NO_DIMS]
Definition: quadtree.hpp:73
QuadTree(double *inp_data, int N)
Definition: quadtree.hpp:98
void fill(int N)
Definition: quadtree.hpp:424
QuadTree & operator=(const QuadTree &)
void init(QuadTree *inp_parent, double *inp_data, double inp_x, double inp_y, double inp_hw, double inp_hh)
Definition: quadtree.hpp:405
void subdivide()
Definition: quadtree.hpp:224
double * data
Definition: quadtree.hpp:85
QuadTree(double *inp_data, double inp_x, double inp_y, double inp_hw, double inp_hh)
Definition: quadtree.hpp:123
double x
Definition: quadtree.hpp:48
void computeNonEdgeForces(int point_index, double theta, double neg_f[], double *sum_Q)
Definition: quadtree.hpp:315