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;
69 static const int QT_NO_DIMS = 2;
70 static const int QT_NODE_CAPACITY = 1;
73 double buff[QT_NO_DIMS];
86 double center_of_mass[QT_NO_DIMS];
87 int index[QT_NODE_CAPACITY];
99 parent(NULL), is_leaf(false), size(0), cum_size(0), boundary(), data(NULL),
100 northWest(NULL), northEast(NULL), southWest(NULL), southEast(NULL)
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];
113 for(
int d = 0; d < QT_NO_DIMS; d++) mean_Y[d] /= (
double) N;
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);
119 delete[] mean_Y;
delete[] max_Y;
delete[] min_Y;
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)
127 init(NULL, inp_data, inp_x, inp_y, inp_hw, inp_hh);
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)
135 init(NULL, inp_data, inp_x, inp_y, inp_hw, inp_hh);
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)
144 init(inp_parent, inp_data, inp_x, inp_y, inp_hw, inp_hh);
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)
153 init(inp_parent, inp_data, inp_x, inp_y, inp_hw, inp_hh);
181 double* point = data + new_index * QT_NO_DIMS;
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];
193 if(is_leaf && size < QT_NODE_CAPACITY) {
194 index[size] = new_index;
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; }
206 any_duplicate = any_duplicate | duplicate;
208 if(any_duplicate)
return true;
211 if(is_leaf) subdivide();
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;
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);
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]);
250 for(
int n = 0; n < size; n++) {
251 double* point = data + index[n] * QT_NO_DIMS;
254 if(!is_leaf)
return northWest->
isCorrect() &&
264 for(
int n = 0; n < size; n++) {
266 double* point = data + index[n] * QT_NO_DIMS;
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;
279 for(
int d = 0; d < QT_NO_DIMS; d++) {
283 if(node->
getParent() == NULL) done =
true;
302 getAllIndices(indices, 0);
307 if(is_leaf)
return 1;
308 return 1 + std::max(std::max(northWest->
getDepth(),
319 if(cum_size == 0 || (is_leaf && size == 1 && index[0] == point_index))
return;
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];
329 if(is_leaf || std::max(boundary.
hh, boundary.
hw)/sqrt(D) < theta) {
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];
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++) {
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);
366 for(
int d = 0; d < QT_NO_DIMS; d++) pos_f[ind1 + d] += D * buff[d];
375 printf(
"Empty node\n");
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");
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");
405 void init(
QuadTree* inp_parent,
double* inp_data,
double inp_x,
double inp_y,
double inp_hw,
double inp_hh)
414 boundary.
hw = inp_hw;
415 boundary.
hh = inp_hh;
420 for(
int i = 0; i < QT_NO_DIMS; i++) center_of_mass[i] = .0;
426 for(
int i = 0; i < N; i++) insert(i);
434 for(
int i = 0; i < size; i++) indices[loc + i] = index[i];
Namespace containing implementation of t-SNE algorithm.
void getAllIndices(int *indices)
bool containsPoint(double point[])
double center_of_mass[QT_NO_DIMS]
void setData(double *inp_data)
void computeEdgeForces(int *row_P, int *col_P, double *val_P, int N, double *pos_f)
QuadTree(double *inp_data, int N, double inp_x, double inp_y, double inp_hw, double inp_hh)
QuadTree(QuadTree *inp_parent, double *inp_data, double inp_x, double inp_y, double inp_hw, double inp_hh)
bool insert(int new_index)
int getAllIndices(int *indices, int loc)
QuadTree(QuadTree *inp_parent, double *inp_data, int N, double inp_x, double inp_y, double inp_hw, double inp_hh)
QuadTree(double *inp_data, int N)
void init(QuadTree *inp_parent, double *inp_data, double inp_x, double inp_y, double inp_hw, double inp_hh)
QuadTree(double *inp_data, double inp_x, double inp_y, double inp_hw, double inp_hh)
void computeNonEdgeForces(int point_index, double theta, double neg_f[], double *sum_Q)