57 if(x - hw > point[0])
return false;
58 if(x + hw < point[0])
return false;
59 if(y - hh > point[1])
return false;
60 if(y + hh < point[1])
return false;
71 static const int QT_NO_DIMS = 2;
72 static const int QT_NODE_CAPACITY = 1;
89 int index[QT_NODE_CAPACITY];
101 parent(NULL), is_leaf(false), size(0), cum_size(0), boundary(), data(NULL),
102 northWest(NULL), northEast(NULL), southWest(NULL), southEast(NULL)
105 ScalarType* mean_Y =
new ScalarType[QT_NO_DIMS];
for(
int d = 0; d < QT_NO_DIMS; d++) mean_Y[d] = .0;
106 ScalarType* min_Y =
new ScalarType[QT_NO_DIMS];
for(
int d = 0; d < QT_NO_DIMS; d++) min_Y[d] = DBL_MAX;
107 ScalarType* max_Y =
new ScalarType[QT_NO_DIMS];
for(
int d = 0; d < QT_NO_DIMS; d++) max_Y[d] = -DBL_MAX;
108 for(
int n = 0; n < N; n++) {
109 for(
int d = 0; d < QT_NO_DIMS; d++) {
110 mean_Y[d] += inp_data[n * QT_NO_DIMS + d];
111 if(inp_data[n * QT_NO_DIMS + d] < min_Y[d]) min_Y[d] = inp_data[n * QT_NO_DIMS + d];
112 if(inp_data[n * QT_NO_DIMS + d] > max_Y[d]) max_Y[d] = inp_data[n * QT_NO_DIMS + d];
115 for(
int d = 0; d < QT_NO_DIMS; d++) mean_Y[d] /= (
ScalarType) N;
118 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,
119 std::max(max_Y[1] - mean_Y[1], mean_Y[1] - min_Y[1]) + 1e-5);
121 delete[] mean_Y;
delete[] max_Y;
delete[] min_Y;
126 parent(NULL), is_leaf(false), size(0), cum_size(0), boundary(), data(NULL),
127 northWest(NULL), northEast(NULL), southWest(NULL), southEast(NULL)
129 init(NULL, inp_data, inp_x, inp_y, inp_hw, inp_hh);
134 parent(NULL), is_leaf(false), size(0), cum_size(0), boundary(), data(NULL),
135 northWest(NULL), northEast(NULL), southWest(NULL), southEast(NULL)
137 init(NULL, inp_data, inp_x, inp_y, inp_hw, inp_hh);
143 parent(NULL), is_leaf(false), size(0), cum_size(0), boundary(), data(NULL),
144 northWest(NULL), northEast(NULL), southWest(NULL), southEast(NULL)
146 init(inp_parent, inp_data, inp_x, inp_y, inp_hw, inp_hh);
152 parent(NULL), is_leaf(false), size(0), cum_size(0), boundary(), data(NULL),
153 northWest(NULL), northEast(NULL), southWest(NULL), southEast(NULL)
155 init(inp_parent, inp_data, inp_x, inp_y, inp_hw, inp_hh);
183 ScalarType* point = data + new_index * QT_NO_DIMS;
191 for(
int d = 0; d < QT_NO_DIMS; d++) center_of_mass[d] *= mult1;
192 for(
int d = 0; d < QT_NO_DIMS; d++) center_of_mass[d] += mult2 * point[d];
195 if(is_leaf && size < QT_NODE_CAPACITY) {
196 index[size] = new_index;
202 bool any_duplicate =
false;
203 for(
int n = 0; n < size; n++) {
204 bool duplicate =
true;
205 for(
int d = 0; d < QT_NO_DIMS; d++) {
206 if(point[d] != data[index[n] * QT_NO_DIMS + d]) { duplicate =
false;
break; }
208 any_duplicate = any_duplicate | duplicate;
210 if(any_duplicate)
return true;
213 if(is_leaf) subdivide();
216 if(northWest->
insert(new_index))
return true;
217 if(northEast->
insert(new_index))
return true;
218 if(southWest->
insert(new_index))
return true;
219 if(southEast->
insert(new_index))
return true;
229 northWest =
new QuadTree(
this, data, boundary.
x - .5 * boundary.
hw, boundary.
y - .5 * boundary.
hh, .5 * boundary.
hw, .5 * boundary.
hh);
230 northEast =
new QuadTree(
this, data, boundary.
x + .5 * boundary.
hw, boundary.
y - .5 * boundary.
hh, .5 * boundary.
hw, .5 * boundary.
hh);
231 southWest =
new QuadTree(
this, data, boundary.
x - .5 * boundary.
hw, boundary.
y + .5 * boundary.
hh, .5 * boundary.
hw, .5 * boundary.
hh);
232 southEast =
new QuadTree(
this, data, boundary.
x + .5 * boundary.
hw, boundary.
y + .5 * boundary.
hh, .5 * boundary.
hw, .5 * boundary.
hh);
235 for(
int i = 0; i < size; i++) {
236 bool success =
false;
237 if(!success) success = northWest->
insert(index[i]);
238 if(!success) success = northEast->
insert(index[i]);
239 if(!success) success = southWest->
insert(index[i]);
240 if(!success) success = southEast->
insert(index[i]);
252 for(
int n = 0; n < size; n++) {
253 ScalarType* point = data + index[n] * QT_NO_DIMS;
256 if(!is_leaf)
return northWest->
isCorrect() &&
266 for(
int n = 0; n < size; n++) {
268 ScalarType* point = data + index[n] * QT_NO_DIMS;
272 int rem_index = index[n];
273 for(
int m = n + 1; m < size; m++) index[m - 1] = index[m];
274 index[size - 1] = -1;
281 for(
int d = 0; d < QT_NO_DIMS; d++) {
285 if(node->
getParent() == NULL) done =
true;
304 getAllIndices(indices, 0);
309 if(is_leaf)
return 1;
310 return 1 + std::max(std::max(northWest->
getDepth(),
321 if(cum_size == 0 || (is_leaf && size == 1 && index[0] == point_index))
return;
325 int ind = point_index * QT_NO_DIMS;
326 for(
int d = 0; d < QT_NO_DIMS; d++) buff[d] = data[ind + d];
327 for(
int d = 0; d < QT_NO_DIMS; d++) buff[d] -= center_of_mass[d];
328 for(
int d = 0; d < QT_NO_DIMS; d++) D += buff[d] * buff[d];
331 if(is_leaf || std::max(boundary.
hh, boundary.
hw)/sqrt(D) < theta) {
335 *sum_Q += cum_size * Q;
337 for(
int d = 0; d < QT_NO_DIMS; d++) neg_f[d] += mult * buff[d];
355 for(
int n = 0; n < N; n++) {
356 ind1 = n * QT_NO_DIMS;
357 for(
int i = row_P[n]; i < row_P[n + 1]; i++) {
361 ind2 = col_P[i] * QT_NO_DIMS;
362 for(
int d = 0; d < QT_NO_DIMS; d++) buff[d] = data[ind1 + d];
363 for(
int d = 0; d < QT_NO_DIMS; d++) buff[d] -= data[ind2 + d];
364 for(
int d = 0; d < QT_NO_DIMS; d++) D += buff[d] * buff[d];
365 D = val_P[i] / (1.0 + D);
368 for(
int d = 0; d < QT_NO_DIMS; d++) pos_f[ind1 + d] += D * buff[d];
377 printf(
"Empty node\n");
382 printf(
"Leaf node; data = [");
383 for(
int i = 0; i < size; i++) {
384 ScalarType* point = data + index[i] * QT_NO_DIMS;
385 for(
int d = 0; d < QT_NO_DIMS; d++) printf(
"%f, ", point[d]);
386 printf(
" (index = %d)", index[i]);
387 if(i < size - 1) printf(
"\n");
392 printf(
"Intersection node with center-of-mass = [");
393 for(
int d = 0; d < QT_NO_DIMS; d++) printf(
"%f, ", center_of_mass[d]);
394 printf(
"]; children are:\n");
416 boundary.
hw = inp_hw;
417 boundary.
hh = inp_hh;
422 for(
int i = 0; i < QT_NO_DIMS; i++) center_of_mass[i] = .0;
428 for(
int i = 0; i < N; i++) insert(i);
436 for(
int i = 0; i < size; i++) indices[loc + i] = index[i];
QuadTree(ScalarType *inp_data, int N)
Namespace containing implementation of t-SNE algorithm.
void getAllIndices(int *indices)
QuadTree(ScalarType *inp_data, int N, ScalarType inp_x, ScalarType inp_y, ScalarType inp_hw, ScalarType inp_hh)
double ScalarType
default scalar value (can be overrided with TAPKEE_CUSTOM_INTERNAL_NUMTYPE define) ...
void computeNonEdgeForces(int point_index, ScalarType theta, ScalarType neg_f[], ScalarType *sum_Q)
QuadTree(QuadTree *inp_parent, ScalarType *inp_data, ScalarType inp_x, ScalarType inp_y, ScalarType inp_hw, ScalarType inp_hh)
QuadTree(QuadTree *inp_parent, ScalarType *inp_data, int N, ScalarType inp_x, ScalarType inp_y, ScalarType inp_hw, ScalarType inp_hh)
void init(QuadTree *inp_parent, ScalarType *inp_data, ScalarType inp_x, ScalarType inp_y, ScalarType inp_hw, ScalarType inp_hh)
QuadTree(ScalarType *inp_data, ScalarType inp_x, ScalarType inp_y, ScalarType inp_hw, ScalarType inp_hh)
void setData(ScalarType *inp_data)
bool insert(int new_index)
int getAllIndices(int *indices, int loc)
void computeEdgeForces(int *row_P, int *col_P, ScalarType *val_P, int N, ScalarType *pos_f)
ScalarType center_of_mass[QT_NO_DIMS]
bool containsPoint(ScalarType point[])