53 static inline double sign(
double x) {
return (x == .0 ? .0 : (x < .0 ? -1.0 : 1.0)); }
58 void run(
double* X,
int N,
int D,
double* Y,
int no_dims,
double perplexity,
double theta)
61 bool exact = (theta == .0) ?
true :
false;
68 float total_time = .0;
70 int max_iter = 1000, stop_lying_iter = 250, mom_switch_iter = 250;
71 double momentum = .5, final_momentum = .8;
75 double* dY = (
double*) malloc(N * no_dims *
sizeof(
double));
76 double* uY = (
double*) malloc(N * no_dims *
sizeof(
double));
77 double* gains = (
double*) malloc(N * no_dims *
sizeof(
double));
78 if(dY == NULL || uY == NULL || gains == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
79 for(
int i = 0; i < N * no_dims; i++) uY[i] = .0;
80 for(
int i = 0; i < N * no_dims; i++) gains[i] = 1.0;
83 double* P=NULL;
int* row_P;
int* col_P;
double* val_P;
89 for(
int i = 0; i < N * D; i++) {
90 if(X[i] > max_X) max_X = X[i];
92 for(
int i = 0; i < N * D; i++) X[i] /= max_X;
98 P = (
double*) malloc(N * N *
sizeof(
double));
99 if(P == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
103 for(
int n = 0; n < N; n++) {
104 for(
int m = n + 1; m < N; m++) {
105 P[n * N + m] += P[m * N + n];
106 P[m * N + n] = P[n * N + m];
110 for(
int i = 0; i < N * N; i++) sum_P += P[i];
111 for(
int i = 0; i < N * N; i++) P[i] /= sum_P;
123 for(
int i = 0; i < row_P[N]; i++) sum_P += val_P[i];
124 for(
int i = 0; i < row_P[N]; i++) val_P[i] /= sum_P;
128 if(exact) {
for(
int i = 0; i < N * N; i++) P[i] *= 12.0; }
129 else {
for(
int i = 0; i < row_P[N]; i++) val_P[i] *= 12.0; }
137 for(
int iter = 0; iter < max_iter; iter++) {
141 else computeGradient(P, row_P, col_P, val_P, Y, N, no_dims, dY, theta);
144 for(
int i = 0; i < N * no_dims; i++) gains[i] = (
sign(dY[i]) !=
sign(uY[i])) ? (gains[i] + .2) : (gains[i] * .8);
145 for(
int i = 0; i < N * no_dims; i++)
if(gains[i] < .01) gains[i] = .01;
148 for(
int i = 0; i < N * no_dims; i++) uY[i] = momentum * uY[i] - eta * gains[i] * dY[i];
149 for(
int i = 0; i < N * no_dims; i++) Y[i] = Y[i] + uY[i];
155 if(iter == stop_lying_iter) {
156 if(exact) {
for(
int i = 0; i < N * N; i++) P[i] /= 12.0; }
157 else {
for(
int i = 0; i < row_P[N]; i++) val_P[i] /= 12.0; }
159 if(iter == mom_switch_iter) momentum = final_momentum;
181 end = clock(); total_time += (float) (end - start) / CLOCKS_PER_SEC;
189 free(row_P); row_P = NULL;
190 free(col_P); col_P = NULL;
191 free(val_P); val_P = NULL;
199 int* row_P = *_row_P;
200 int* col_P = *_col_P;
201 double* val_P = *_val_P;
204 int* row_counts = (
int*) calloc(N,
sizeof(
int));
205 if(row_counts == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
206 for(
int n = 0; n < N; n++) {
207 for(
int i = row_P[n]; i < row_P[n + 1]; i++) {
209 bool present =
false;
210 for(
int m = row_P[col_P[i]]; m < row_P[col_P[i] + 1]; m++) {
211 if(col_P[m] == n) present =
true;
213 if(present) row_counts[n]++;
216 row_counts[col_P[i]]++;
221 for(
int n = 0; n < N; n++) no_elem += row_counts[n];
224 int* sym_row_P = (
int*) malloc((N + 1) *
sizeof(int));
225 int* sym_col_P = (
int*) malloc(no_elem *
sizeof(
int));
226 double* sym_val_P = (
double*) malloc(no_elem *
sizeof(
double));
227 if(sym_row_P == NULL || sym_col_P == NULL || sym_val_P == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
231 for(
int n = 0; n < N; n++) sym_row_P[n + 1] = sym_row_P[n] + row_counts[n];
234 int* offset = (
int*) calloc(N,
sizeof(
int));
235 if(offset == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
236 for(
int n = 0; n < N; n++) {
237 for(
int i = row_P[n]; i < row_P[n + 1]; i++) {
240 bool present =
false;
241 for(
int m = row_P[col_P[i]]; m < row_P[col_P[i] + 1]; m++) {
245 sym_col_P[sym_row_P[n] + offset[n]] = col_P[i];
246 sym_col_P[sym_row_P[col_P[i]] + offset[col_P[i]]] = n;
247 sym_val_P[sym_row_P[n] + offset[n]] = val_P[i] + val_P[m];
248 sym_val_P[sym_row_P[col_P[i]] + offset[col_P[i]]] = val_P[i] + val_P[m];
255 sym_col_P[sym_row_P[n] + offset[n]] = col_P[i];
256 sym_col_P[sym_row_P[col_P[i]] + offset[col_P[i]]] = n;
257 sym_val_P[sym_row_P[n] + offset[n]] = val_P[i];
258 sym_val_P[sym_row_P[col_P[i]] + offset[col_P[i]]] = val_P[i];
262 if(!present || (present && n <= col_P[i])) {
264 if(col_P[i] != n) offset[col_P[i]]++;
270 for(
int i = 0; i < no_elem; i++) sym_val_P[i] /= 2.0;
273 free(*_row_P); *_row_P = sym_row_P;
274 free(*_col_P); *_col_P = sym_col_P;
275 free(*_val_P); *_val_P = sym_val_P;
278 free(offset); offset = NULL;
279 free(row_counts); row_counts = NULL;
284 void computeGradient(
double* ,
int* inp_row_P,
int* inp_col_P,
double* inp_val_P,
double* Y,
int N,
int D,
double* dC,
double theta)
291 double* pos_f = (
double*) calloc(N * D,
sizeof(
double));
292 double* neg_f = (
double*) calloc(N * D,
sizeof(
double));
293 if(pos_f == NULL || neg_f == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
298 for(
int i = 0; i < N * D; i++) {
299 dC[i] = pos_f[i] - (neg_f[i] / sum_Q);
309 for(
int i = 0; i < N * D; i++) dC[i] = 0.0;
312 double* DD = (
double*) malloc(N * N *
sizeof(
double));
313 if(DD == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
317 double* Q = (
double*) malloc(N * N *
sizeof(
double));
318 if(Q == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
320 for(
int n = 0; n < N; n++) {
321 for(
int m = 0; m < N; m++) {
323 Q[n * N + m] = 1 / (1 + DD[n * N + m]);
324 sum_Q += Q[n * N + m];
330 for(
int n = 0; n < N; n++) {
331 for(
int m = 0; m < N; m++) {
333 double mult = (P[n * N + m] - (Q[n * N + m] / sum_Q)) * Q[n * N + m];
334 for(
int d = 0; d < D; d++) {
335 dC[n * D + d] += (Y[n * D + d] - Y[m * D + d]) * mult;
349 double* DD = (
double*) malloc(N * N *
sizeof(
double));
350 double* Q = (
double*) malloc(N * N *
sizeof(
double));
351 if(DD == NULL || Q == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
355 double sum_Q = DBL_MIN;
356 for(
int n = 0; n < N; n++) {
357 for(
int m = 0; m < N; m++) {
359 Q[n * N + m] = 1 / (1 + DD[n * N + m]);
360 sum_Q += Q[n * N + m];
362 else Q[n * N + m] = DBL_MIN;
365 for(
int i = 0; i < N * N; i++) Q[i] /= sum_Q;
369 for(
int n = 0; n < N; n++) {
370 for(
int m = 0; m < N; m++) {
371 C += P[n * N + m] * log((P[n * N + m] + 1e-9) / (Q[n * N + m] + 1e-9));
381 double evaluateError(
int* row_P,
int* col_P,
double* val_P,
double* Y,
int N,
double theta)
384 const int QT_NO_DIMS = 2;
386 double buff[QT_NO_DIMS] = {.0, .0};
393 for(
int n = 0; n < N; n++) {
394 ind1 = n * QT_NO_DIMS;
395 for(
int i = row_P[n]; i < row_P[n + 1]; i++) {
397 ind2 = col_P[i] * QT_NO_DIMS;
398 for(
int d = 0; d < QT_NO_DIMS; d++) buff[d] = Y[ind1 + d];
399 for(
int d = 0; d < QT_NO_DIMS; d++) buff[d] -= Y[ind2 + d];
400 for(
int d = 0; d < QT_NO_DIMS; d++) Q += buff[d] * buff[d];
401 Q = (1.0 / (1.0 + Q)) / sum_Q;
402 C += val_P[i] * log((val_P[i] + FLT_MIN) / (Q + FLT_MIN));
411 double* mean = (
double*) calloc(D,
sizeof(
double));
412 if(mean == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
413 for(
int n = 0; n < N; n++) {
414 for(
int d = 0; d < D; d++) {
415 mean[d] += X[n * D + d];
418 for(
int d = 0; d < D; d++) {
419 mean[d] /= (double) N;
423 for(
int n = 0; n < N; n++) {
424 for(
int d = 0; d < D; d++) {
425 X[n * D + d] -= mean[d];
428 free(mean); mean = NULL;
434 double* DD = (
double*) malloc(N * N *
sizeof(
double));
435 if(DD == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
439 for(
int n = 0; n < N; n++) {
444 double min_beta = -DBL_MAX;
445 double max_beta = DBL_MAX;
451 while(!found && iter < 200) {
454 for(
int m = 0; m < N; m++) P[n * N + m] = exp(-beta * DD[n * N + m]);
455 P[n * N + n] = DBL_MIN;
459 for(
int m = 0; m < N; m++) sum_P += P[n * N + m];
461 for(
int m = 0; m < N; m++) H += beta * (DD[n * N + m] * P[n * N + m]);
462 H = (H / sum_P) + log(sum_P);
465 double Hdiff = H - log(perplexity);
466 if(Hdiff < tol && -Hdiff < tol) {
472 if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
475 beta = (beta + max_beta) / 2.0;
479 if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
482 beta = (beta + min_beta) / 2.0;
491 for(
int m = 0; m < N; m++) P[n * N + m] /= sum_P;
500 if(perplexity > K) printf(
"Perplexity should be lower than K!\n");
503 *_row_P = (
int*) malloc((N + 1) *
sizeof(int));
504 *_col_P = (
int*) calloc(N * K,
sizeof(
int));
505 *_val_P = (
double*) calloc(N * K,
sizeof(
double));
506 if(*_row_P == NULL || *_col_P == NULL || *_val_P == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
507 int* row_P = *_row_P;
508 int* col_P = *_col_P;
509 double* val_P = *_val_P;
510 double* cur_P = (
double*) malloc((N - 1) *
sizeof(double));
511 if(cur_P == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
513 for(
int n = 0; n < N; n++) row_P[n + 1] = row_P[n] + K;
517 std::vector<DataPoint> obj_X(N,
DataPoint(D, -1, X));
518 for(
int n = 0; n < N; n++) obj_X[n] =
DataPoint(D, n, X + n * D);
523 std::vector<DataPoint> indices;
524 std::vector<double> distances;
525 for(
int n = 0; n < N; n++) {
532 tree->
search(obj_X[n], K + 1, &indices, &distances);
537 double min_beta = -DBL_MAX;
538 double max_beta = DBL_MAX;
542 int iter = 0;
double sum_P;
543 while(!found && iter < 200) {
546 for(
int m = 0; m < K; m++) cur_P[m] = exp(-beta * distances[m + 1]);
550 for(
int m = 0; m < K; m++) sum_P += cur_P[m];
552 for(
int m = 0; m < K; m++) H += beta * (distances[m + 1] * cur_P[m]);
553 H = (H / sum_P) + log(sum_P);
556 double Hdiff = H - log(perplexity);
557 if(Hdiff < tol && -Hdiff < tol) {
563 if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
566 beta = (beta + max_beta) / 2.0;
570 if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
573 beta = (beta + min_beta) / 2.0;
582 for(
int m = 0; m < K; m++) cur_P[m] /= sum_P;
583 for(
int m = 0; m < K; m++) {
584 col_P[row_P[n] + m] = indices[m + 1].index();
585 val_P[row_P[n] + m] = cur_P[m];
595 void computeGaussianPerplexity(
double* X,
int N,
int D,
int** _row_P,
int** _col_P,
double** _val_P,
double perplexity,
double threshold)
598 double* buff = (
double*) malloc(D *
sizeof(
double));
599 double* DD = (
double*) malloc(N *
sizeof(
double));
600 double* cur_P = (
double*) malloc(N *
sizeof(
double));
601 if(buff == NULL || DD == NULL || cur_P == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
605 for(
int n = 0; n < N; n++) {
608 for(
int m = 0; m < N; m++) {
609 for(
int d = 0; d < D; d++) buff[d] = X[n * D + d];
610 for(
int d = 0; d < D; d++) buff[d] -= X[m * D + d];
612 for(
int d = 0; d < D; d++) DD[m] += buff[d] * buff[d];
618 double min_beta = -DBL_MAX;
619 double max_beta = DBL_MAX;
623 int iter = 0;
double sum_P;
624 while(!found && iter < 200) {
627 for(
int m = 0; m < N; m++) cur_P[m] = exp(-beta * DD[m]);
632 for(
int m = 0; m < N; m++) sum_P += cur_P[m];
634 for(
int m = 0; m < N; m++) H += beta * (DD[m] * cur_P[m]);
635 H = (H / sum_P) + log(sum_P);
638 double Hdiff = H - log(perplexity);
639 if(Hdiff < tol && -Hdiff < tol) {
645 if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
648 beta = (beta + max_beta) / 2.0;
652 if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
655 beta = (beta + min_beta) / 2.0;
664 for(
int m = 0; m < N; m++) cur_P[m] /= sum_P;
665 for(
int m = 0; m < N; m++) {
666 if(cur_P[m] > threshold / (
double) N) total_count++;
671 *_row_P = (
int*) malloc((N + 1) *
sizeof(int));
672 *_col_P = (
int*) malloc(total_count *
sizeof(
int));
673 *_val_P = (
double*) malloc(total_count *
sizeof(
double));
674 int* row_P = *_row_P;
675 int* col_P = *_col_P;
676 double* val_P = *_val_P;
677 if(row_P == NULL || col_P == NULL || val_P == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
682 for(
int n = 0; n < N; n++) {
685 for(
int m = 0; m < N; m++) {
686 for(
int d = 0; d < D; d++) buff[d] = X[n * D + d];
687 for(
int d = 0; d < D; d++) buff[d] -= X[m * D + d];
689 for(
int d = 0; d < D; d++) DD[m] += buff[d] * buff[d];
695 double min_beta = -DBL_MAX;
696 double max_beta = DBL_MAX;
700 int iter = 0;
double sum_P;
701 while(!found && iter < 200) {
704 for(
int m = 0; m < N; m++) cur_P[m] = exp(-beta * DD[m]);
709 for(
int m = 0; m < N; m++) sum_P += cur_P[m];
711 for(
int m = 0; m < N; m++) H += beta * (DD[m] * cur_P[m]);
712 H = (H / sum_P) + log(sum_P);
715 double Hdiff = H - log(perplexity);
716 if(Hdiff < tol && -Hdiff < tol) {
722 if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
725 beta = (beta + max_beta) / 2.0;
729 if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
732 beta = (beta + min_beta) / 2.0;
741 for(
int m = 0; m < N; m++) cur_P[m] /= sum_P;
742 for(
int m = 0; m < N; m++) {
743 if(cur_P[m] > threshold / (
double) N) {
745 val_P[count] = cur_P[m];
749 row_P[n + 1] = count;
754 free(buff); buff = NULL;
755 free(cur_P); cur_P = NULL;
760 double* dataSums = (
double*) calloc(N,
sizeof(
double));
761 if(dataSums == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
762 for(
int n = 0; n < N; n++) {
763 for(
int d = 0; d < D; d++) {
764 dataSums[n] += (X[n * D + d] * X[n * D + d]);
767 for(
int n = 0; n < N; n++) {
768 for(
int m = 0; m < N; m++) {
769 DD[n * N + m] = dataSums[n] + dataSums[m];
772 Eigen::Map<Eigen::MatrixXd> DD_map(DD,N,N);
773 Eigen::Map<Eigen::MatrixXd> X_map(X,D,N);
774 DD_map.noalias() = -2.0*X_map.transpose()*X_map;
777 free(dataSums); dataSums = NULL;
void run(double *X, int N, int D, double *Y, int no_dims, double perplexity, double theta)
void symmetrizeMatrix(int **_row_P, int **_col_P, double **_val_P, int N)
void search(const T &target, int k, std::vector< T > *results, std::vector< double > *distances)
void computeGaussianPerplexity(double *X, int N, int D, int **_row_P, int **_col_P, double **_val_P, double perplexity, int K)
double evaluateError(double *P, double *Y, int N)
void computeGradient(double *, int *inp_row_P, int *inp_col_P, double *inp_val_P, double *Y, int N, int D, double *dC, double theta)
double evaluateError(int *row_P, int *col_P, double *val_P, double *Y, int N, double theta)
void computeEdgeForces(int *row_P, int *col_P, double *val_P, int N, double *pos_f)
void computeGaussianPerplexity(double *X, int N, int D, int **_row_P, int **_col_P, double **_val_P, double perplexity, double threshold)
static double sign(double x)
void create(const std::vector< T > &items)
void computeGaussianPerplexity(double *X, int N, int D, double *P, double perplexity)
void message_info(const std::string &msg)
static LoggingSingleton & instance()
ScalarType gaussian_random()
void computeExactGradient(double *P, double *Y, int N, int D, double *dC)
void computeSquaredEuclideanDistance(double *X, int N, int D, double *DD)
void zeroMean(double *X, int N, int D)
void computeNonEdgeForces(int point_index, double theta, double neg_f[], double *sum_Q)