63 bool exact = (theta == .0) ?
true :
false;
70 float total_time = .0;
71 int max_iter = 1000, stop_lying_iter = 250, mom_switch_iter = 250;
80 gains.setConstant(1.0);
83 int* row_P=NULL;
int* col_P=NULL;
ScalarType* val_P=NULL;
87 X.array() /= X.maxCoeff();
97 for(
int n = 0; n < N; n++) {
98 for(
int m = n + 1; m < N; m++) {
99 P.data()[n * N + m] += P.data()[m * N + n];
100 P.data()[m * N + n] = P.data()[n * N + m];
103 P.array() /= P.array().sum();
115 for(
int i = 0; i < row_P[N]; i++) sum_P += val_P[i];
116 for(
int i = 0; i < row_P[N]; i++) val_P[i] /= sum_P;
124 for(
int i = 0; i < row_P[N]; i++) val_P[i] *= 12.0;
133 for(
int iter = 0; iter < max_iter; iter++) {
137 else computeGradient(P.data(), row_P, col_P, val_P, Y, N, no_dims, dY.data(), theta);
140 for(
int i = 0; i < N * no_dims; i++) gains.data()[i] = (
sign(dY.data()[i]) !=
sign(uY.data()[i])) ? (gains.data()[i] + .2) : (gains.data()[i] * .8);
141 for(
int i = 0; i < N * no_dims; i++)
if(gains.data()[i] < .01) gains.data()[i] = .01;
144 for(
int i = 0; i < N * no_dims; i++) uY.data()[i] = momentum * uY.data()[i] - eta * gains.data()[i] * dY.data()[i];
145 for(
int i = 0; i < N * no_dims; i++) Y[i] = Y[i] + uY.data()[i];
151 if(iter == stop_lying_iter) {
156 for(
int i = 0; i < row_P[N]; i++) val_P[i] /= 12.0;
159 if(iter == mom_switch_iter) momentum = final_momentum;
162 if((iter > 0) && ((iter % 50 == 0) || (iter == max_iter - 1))) {
173 free(row_P); row_P = NULL;
174 free(col_P); col_P = NULL;
175 free(val_P); val_P = NULL;
183 int* row_P = *_row_P;
184 int* col_P = *_col_P;
188 int* row_counts = (
int*) calloc(N,
sizeof(
int));
189 if(row_counts == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
190 for(
int n = 0; n < N; n++) {
191 for(
int i = row_P[n]; i < row_P[n + 1]; i++) {
193 bool present =
false;
194 for(
int m = row_P[col_P[i]]; m < row_P[col_P[i] + 1]; m++) {
195 if(col_P[m] == n) present =
true;
197 if(present) row_counts[n]++;
200 row_counts[col_P[i]]++;
205 for(
int n = 0; n < N; n++) no_elem += row_counts[n];
208 int* sym_row_P = (
int*) malloc((N + 1) *
sizeof(int));
209 int* sym_col_P = (
int*) malloc(no_elem *
sizeof(
int));
211 if(sym_row_P == NULL || sym_col_P == NULL || sym_val_P == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
215 for(
int n = 0; n < N; n++) sym_row_P[n + 1] = sym_row_P[n] + row_counts[n];
218 int* offset = (
int*) calloc(N,
sizeof(
int));
219 if(offset == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
220 for(
int n = 0; n < N; n++) {
221 for(
int i = row_P[n]; i < row_P[n + 1]; i++) {
224 bool present =
false;
225 for(
int m = row_P[col_P[i]]; m < row_P[col_P[i] + 1]; m++) {
229 sym_col_P[sym_row_P[n] + offset[n]] = col_P[i];
230 sym_col_P[sym_row_P[col_P[i]] + offset[col_P[i]]] = n;
231 sym_val_P[sym_row_P[n] + offset[n]] = val_P[i] + val_P[m];
232 sym_val_P[sym_row_P[col_P[i]] + offset[col_P[i]]] = val_P[i] + val_P[m];
239 sym_col_P[sym_row_P[n] + offset[n]] = col_P[i];
240 sym_col_P[sym_row_P[col_P[i]] + offset[col_P[i]]] = n;
241 sym_val_P[sym_row_P[n] + offset[n]] = val_P[i];
242 sym_val_P[sym_row_P[col_P[i]] + offset[col_P[i]]] = val_P[i];
246 if(!present || (present && n <= col_P[i])) {
248 if(col_P[i] != n) offset[col_P[i]]++;
254 for(
int i = 0; i < no_elem; i++) sym_val_P[i] /= 2.0;
257 free(*_row_P); *_row_P = sym_row_P;
258 free(*_col_P); *_col_P = sym_col_P;
259 free(*_val_P); *_val_P = sym_val_P;
262 free(offset); offset = NULL;
263 free(row_counts); row_counts = NULL;
277 if(pos_f == NULL || neg_f == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
282 for(
int i = 0; i < N * D; i++) {
283 dC[i] = pos_f[i] - (neg_f[i] / sum_Q);
293 for(
int i = 0; i < N * D; i++) dC[i] = 0.0;
297 if(DD == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
302 if(Q == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
304 for(
int n = 0; n < N; n++) {
305 for(
int m = 0; m < N; m++) {
307 Q[n * N + m] = 1 / (1 + DD[n * N + m]);
308 sum_Q += Q[n * N + m];
314 for(
int n = 0; n < N; n++) {
315 for(
int m = 0; m < N; m++) {
317 ScalarType mult = (P[n * N + m] - (Q[n * N + m] / sum_Q)) * Q[n * N + m];
318 for(
int d = 0; d < D; d++) {
319 dC[n * D + d] += (Y[n * D + d] - Y[m * D + d]) * mult;
335 if(DD == NULL || Q == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
340 for(
int n = 0; n < N; n++) {
341 for(
int m = 0; m < N; m++) {
343 Q[n * N + m] = 1 / (1 + DD[n * N + m]);
344 sum_Q += Q[n * N + m];
346 else Q[n * N + m] = DBL_MIN;
349 for(
int i = 0; i < N * N; i++) Q[i] /= sum_Q;
353 for(
int n = 0; n < N; n++) {
354 for(
int m = 0; m < N; m++) {
355 C += P[n * N + m] * log((P[n * N + m] + 1e-9) / (Q[n * N + m] + 1e-9));
368 const int QT_NO_DIMS = 2;
377 for(
int n = 0; n < N; n++) {
378 ind1 = n * QT_NO_DIMS;
379 for(
int i = row_P[n]; i < row_P[n + 1]; i++) {
381 ind2 = col_P[i] * QT_NO_DIMS;
382 for(
int d = 0; d < QT_NO_DIMS; d++) buff[d] = Y[ind1 + d];
383 for(
int d = 0; d < QT_NO_DIMS; d++) buff[d] -= Y[ind2 + d];
384 for(
int d = 0; d < QT_NO_DIMS; d++) Q += buff[d] * buff[d];
385 Q = (1.0 / (1.0 + Q)) / sum_Q;
386 C += val_P[i] * log((val_P[i] + FLT_MIN) / (Q + FLT_MIN));
396 if(mean == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
397 for(
int n = 0; n < N; n++) {
398 for(
int d = 0; d < D; d++) {
399 mean[d] += X[n * D + d];
402 for(
int d = 0; d < D; d++) {
407 for(
int n = 0; n < N; n++) {
408 for(
int d = 0; d < D; d++) {
409 X[n * D + d] -= mean[d];
412 free(mean); mean = NULL;
419 if(DD == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
423 for(
int n = 0; n < N; n++) {
435 while(!found && iter < 200) {
438 for(
int m = 0; m < N; m++) P[n * N + m] = exp(-beta * DD[n * N + m]);
439 P[n * N + n] = DBL_MIN;
443 for(
int m = 0; m < N; m++) sum_P += P[n * N + m];
445 for(
int m = 0; m < N; m++) H += beta * (DD[n * N + m] * P[n * N + m]);
446 H = (H / sum_P) + log(sum_P);
450 if(Hdiff < tol && -Hdiff < tol) {
456 if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
459 beta = (beta + max_beta) / 2.0;
463 if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
466 beta = (beta + min_beta) / 2.0;
475 for(
int m = 0; m < N; m++) P[n * N + m] /= sum_P;
484 if(perplexity > K) printf(
"Perplexity should be lower than K!\n");
487 *_row_P = (
int*) malloc((N + 1) *
sizeof(int));
488 *_col_P = (
int*) calloc(N * K,
sizeof(
int));
490 if(*_row_P == NULL || *_col_P == NULL || *_val_P == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
491 int* row_P = *_row_P;
492 int* col_P = *_col_P;
495 if(cur_P == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
497 for(
int n = 0; n < N; n++) row_P[n + 1] = row_P[n] + K;
501 std::vector<DataPoint> obj_X(N,
DataPoint(D, -1, X));
502 for(
int n = 0; n < N; n++) obj_X[n] =
DataPoint(D, n, X + n * D);
507 std::vector<DataPoint> indices;
508 std::vector<ScalarType> distances;
509 for(
int n = 0; n < N; n++) {
516 tree->
search(obj_X[n], K + 1, &indices, &distances);
527 while(!found && iter < 200) {
530 for(
int m = 0; m < K; m++) cur_P[m] = exp(-beta * distances[m + 1]);
534 for(
int m = 0; m < K; m++) sum_P += cur_P[m];
536 for(
int m = 0; m < K; m++) H += beta * (distances[m + 1] * cur_P[m]);
537 H = (H / sum_P) + log(sum_P);
541 if(Hdiff < tol && -Hdiff < tol) {
547 if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
550 beta = (beta + max_beta) / 2.0;
554 if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
557 beta = (beta + min_beta) / 2.0;
566 for(
int m = 0; m < K; m++) cur_P[m] /= sum_P;
567 for(
int m = 0; m < K; m++) {
568 col_P[row_P[n] + m] = indices[m + 1].index();
569 val_P[row_P[n] + m] = cur_P[m];
585 if(buff == NULL || DD == NULL || cur_P == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
589 for(
int n = 0; n < N; n++) {
592 for(
int m = 0; m < N; m++) {
593 for(
int d = 0; d < D; d++) buff[d] = X[n * D + d];
594 for(
int d = 0; d < D; d++) buff[d] -= X[m * D + d];
596 for(
int d = 0; d < D; d++) DD[m] += buff[d] * buff[d];
608 while(!found && iter < 200) {
611 for(
int m = 0; m < N; m++) cur_P[m] = exp(-beta * DD[m]);
616 for(
int m = 0; m < N; m++) sum_P += cur_P[m];
618 for(
int m = 0; m < N; m++) H += beta * (DD[m] * cur_P[m]);
619 H = (H / sum_P) + log(sum_P);
623 if(Hdiff < tol && -Hdiff < tol) {
629 if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
632 beta = (beta + max_beta) / 2.0;
636 if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
639 beta = (beta + min_beta) / 2.0;
648 for(
int m = 0; m < N; m++) cur_P[m] /= sum_P;
649 for(
int m = 0; m < N; m++) {
650 if(cur_P[m] > threshold / (
ScalarType) N) total_count++;
655 *_row_P = (
int*) malloc((N + 1) *
sizeof(int));
656 *_col_P = (
int*) malloc(total_count *
sizeof(
int));
658 int* row_P = *_row_P;
659 int* col_P = *_col_P;
661 if(row_P == NULL || col_P == NULL || val_P == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
666 for(
int n = 0; n < N; n++) {
669 for(
int m = 0; m < N; m++) {
670 for(
int d = 0; d < D; d++) buff[d] = X[n * D + d];
671 for(
int d = 0; d < D; d++) buff[d] -= X[m * D + d];
673 for(
int d = 0; d < D; d++) DD[m] += buff[d] * buff[d];
685 while(!found && iter < 200) {
688 for(
int m = 0; m < N; m++) cur_P[m] = exp(-beta * DD[m]);
693 for(
int m = 0; m < N; m++) sum_P += cur_P[m];
695 for(
int m = 0; m < N; m++) H += beta * (DD[m] * cur_P[m]);
696 H = (H / sum_P) + log(sum_P);
700 if(Hdiff < tol && -Hdiff < tol) {
706 if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
709 beta = (beta + max_beta) / 2.0;
713 if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
716 beta = (beta + min_beta) / 2.0;
725 for(
int m = 0; m < N; m++) cur_P[m] /= sum_P;
726 for(
int m = 0; m < N; m++) {
729 val_P[count] = cur_P[m];
733 row_P[n + 1] = count;
738 free(buff); buff = NULL;
739 free(cur_P); cur_P = NULL;
745 if(dataSums == NULL) { printf(
"Memory allocation failed!\n"); exit(1); }
746 for(
int n = 0; n < N; n++) {
747 for(
int d = 0; d < D; d++) {
748 dataSums[n] += (X[n * D + d] * X[n * D + d]);
751 for(
int n = 0; n < N; n++) {
752 for(
int m = 0; m < N; m++) {
753 DD[n * N + m] = dataSums[n] + dataSums[m];
756 Eigen::Map<tapkee::DenseMatrix> DD_map(DD,N,N);
757 Eigen::Map<tapkee::DenseMatrix> X_map(X,D,N);
758 DD_map.noalias() = -2.0*X_map.transpose()*X_map;
761 free(dataSums); dataSums = NULL;
Eigen::Matrix< tapkee::ScalarType, Eigen::Dynamic, Eigen::Dynamic > DenseMatrix
dense matrix type (non-overridable)
ScalarType evaluateError(ScalarType *P, ScalarType *Y, int N)
Namespace containing implementation of t-SNE algorithm.
void computeExactGradient(ScalarType *P, ScalarType *Y, int N, int D, ScalarType *dC)
void search(const T &target, int k, std::vector< T > *results, std::vector< ScalarType > *distances)
double ScalarType
default scalar value (can be overrided with TAPKEE_CUSTOM_INTERNAL_NUMTYPE define) ...
void run(tapkee::DenseMatrix &X, int N, int D, ScalarType *Y, int no_dims, ScalarType perplexity, ScalarType theta)
void computeNonEdgeForces(int point_index, ScalarType theta, ScalarType neg_f[], ScalarType *sum_Q)
void computeGaussianPerplexity(ScalarType *X, int N, int D, int **_row_P, int **_col_P, ScalarType **_val_P, ScalarType perplexity, int K)
void computeGradient(ScalarType *, int *inp_row_P, int *inp_col_P, ScalarType *inp_val_P, ScalarType *Y, int N, int D, ScalarType *dC, ScalarType theta)
void zeroMean(ScalarType *X, int N, int D)
void symmetrizeMatrix(int **_row_P, int **_col_P, ScalarType **_val_P, int N)
void computeSquaredEuclideanDistance(ScalarType *X, int N, int D, ScalarType *DD)
void create(const std::vector< T > &items)
void computeEdgeForces(int *row_P, int *col_P, ScalarType *val_P, int N, ScalarType *pos_f)
void message_info(const std::string &msg)
static LoggingSingleton & instance()
ScalarType gaussian_random()
void computeGaussianPerplexity(ScalarType *X, int N, int D, int **_row_P, int **_col_P, ScalarType **_val_P, ScalarType perplexity, ScalarType threshold)
void computeGaussianPerplexity(ScalarType *X, int N, int D, ScalarType *P, ScalarType perplexity)
ScalarType evaluateError(int *row_P, int *col_P, ScalarType *val_P, ScalarType *Y, int N, ScalarType theta)
static ScalarType sign(ScalarType x)