Tapkee
tsne.hpp
Go to the documentation of this file.
1 
32 #ifndef TSNE_H
33 #define TSNE_H
34 
35 /* Tapkee includes */
36 #include <tapkee/utils/logging.hpp>
37 #include <tapkee/utils/time.hpp>
40 /* End of Tapkee includes */
41 
42 #include <math.h>
43 #include <float.h>
44 #include <stdlib.h>
45 #include <stdio.h>
46 #include <cstring>
47 #include <time.h>
48 
50 namespace tsne
51 {
52 
53 using tapkee::ScalarType;
54 
55 static inline ScalarType sign(ScalarType x) { return (x == .0 ? .0 : (x < .0 ? -1.0 : 1.0)); }
56 
57 class TSNE
58 {
59 public:
60  void run(tapkee::DenseMatrix& X, int N, int D, ScalarType* Y, int no_dims, ScalarType perplexity, ScalarType theta)
61  {
62  // Determine whether we are using an exact algorithm
63  bool exact = (theta == .0) ? true : false;
64  if (exact)
65  tapkee::LoggingSingleton::instance().message_info("Using exact t-SNE algorithm");
66  else
67  tapkee::LoggingSingleton::instance().message_info("Using Barnes-Hut-SNE algorithm");
68 
69  // Set learning parameters
70  float total_time = .0;
71  int max_iter = 1000, stop_lying_iter = 250, mom_switch_iter = 250;
72  ScalarType momentum = .5, final_momentum = .8;
73  ScalarType eta = 200.0;
74 
75  tapkee::DenseMatrix dY(N, no_dims);
76  tapkee::DenseMatrix uY(N, no_dims);
77  tapkee::DenseMatrix gains(N, no_dims);
79  uY.setZero();
80  gains.setConstant(1.0);
81 
82  // Normalize input data (to prevent numerical problems)
83  int* row_P=NULL; int* col_P=NULL; ScalarType* val_P=NULL;
84  {
85  tapkee::tapkee_internal::timed_context context("Input similarities computation");
86  zeroMean(X.data(), N, D);
87  X.array() /= X.maxCoeff();
88 
89  // Compute input similarities for exact t-SNE
90  if(exact) {
91 
92  // Compute similarities
93  P = tapkee::DenseMatrix(N,N);
94  computeGaussianPerplexity(X.data(), N, D, P.data(), perplexity);
95 
96  // Symmetrize input similarities
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];
101  }
102  }
103  P.array() /= P.array().sum();
104  }
105 
106  // Compute input similarities for approximate t-SNE
107  else {
108 
109  // Compute asymmetric pairwise input similarities
110  computeGaussianPerplexity(X.data(), N, D, &row_P, &col_P, &val_P, perplexity, (int) (3 * perplexity));
111 
112  // Symmetrize input similarities
113  symmetrizeMatrix(&row_P, &col_P, &val_P, N);
114  ScalarType sum_P = .0;
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;
117  }
118 
119  // Lie about the P-values
120  if(exact)
121  P.array() *= 12.0;
122  else
123  {
124  for(int i = 0; i < row_P[N]; i++) val_P[i] *= 12.0;
125  }
126 
127  // Initialize solution (randomly)
128  for(int i = 0; i < N * no_dims; i++) Y[i] = tapkee::gaussian_random() * .0001;
129  }
130 
131  {
132  tapkee::tapkee_internal::timed_context context("Main t-SNE loop");
133  for(int iter = 0; iter < max_iter; iter++) {
134 
135  // Compute (approximate) gradient
136  if(exact) computeExactGradient(P.data(), Y, N, no_dims, dY.data());
137  else computeGradient(P.data(), row_P, col_P, val_P, Y, N, no_dims, dY.data(), theta);
138 
139  // Update gains
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;
142 
143  // Perform gradient update (with momentum and gains)
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];
146 
147  // Make solution zero-mean
148  zeroMean(Y, N, no_dims);
149 
150  // Stop lying about the P-values after a while, and switch momentum
151  if(iter == stop_lying_iter) {
152  if (exact)
153  P.array() /= 12.0;
154  else
155  {
156  for(int i = 0; i < row_P[N]; i++) val_P[i] /= 12.0;
157  }
158  }
159  if(iter == mom_switch_iter) momentum = final_momentum;
160 
161  // Print out progress
162  if((iter > 0) && ((iter % 50 == 0) || (iter == max_iter - 1))) {
163  ScalarType C = .0;
164  if(exact) C = evaluateError(P.data(), Y, N);
165  else C = evaluateError(row_P, col_P, val_P, Y, N, theta); // doing approximate computation here!
167  formatting::format("Iteration {}: error is {}\n", iter, C));
168  }
169  }
170  // Clean up memory
171  if(!exact)
172  {
173  free(row_P); row_P = NULL;
174  free(col_P); col_P = NULL;
175  free(val_P); val_P = NULL;
176  }
177  }
178  }
179 
180  void symmetrizeMatrix(int** _row_P, int** _col_P, ScalarType** _val_P, int N)
181  {
182  // Get sparse matrix
183  int* row_P = *_row_P;
184  int* col_P = *_col_P;
185  ScalarType* val_P = *_val_P;
186 
187  // Count number of elements and row counts of symmetric matrix
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++) {
192  // Check whether element (col_P[i], n) is present
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;
196  }
197  if(present) row_counts[n]++;
198  else {
199  row_counts[n]++;
200  row_counts[col_P[i]]++;
201  }
202  }
203  }
204  int no_elem = 0;
205  for(int n = 0; n < N; n++) no_elem += row_counts[n];
206 
207  // Allocate memory for symmetrized matrix
208  int* sym_row_P = (int*) malloc((N + 1) * sizeof(int));
209  int* sym_col_P = (int*) malloc(no_elem * sizeof(int));
210  ScalarType* sym_val_P = (ScalarType*) malloc(no_elem * sizeof(ScalarType));
211  if(sym_row_P == NULL || sym_col_P == NULL || sym_val_P == NULL) { printf("Memory allocation failed!\n"); exit(1); }
212 
213  // Construct new row indices for symmetric matrix
214  sym_row_P[0] = 0;
215  for(int n = 0; n < N; n++) sym_row_P[n + 1] = sym_row_P[n] + row_counts[n];
216 
217  // Fill the result matrix
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++) { // considering element(n, col_P[i])
222 
223  // Check whether element (col_P[i], n) is present
224  bool present = false;
225  for(int m = row_P[col_P[i]]; m < row_P[col_P[i] + 1]; m++) {
226  if(col_P[m] == n) {
227  present = true;
228  if(n <= col_P[i]) { // make sure we do not add elements twice
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];
233  }
234  }
235  }
236 
237  // If (col_P[i], n) is not present, there is no addition involved
238  if(!present) {
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];
243  }
244 
245  // Update offsets
246  if(!present || (present && n <= col_P[i])) {
247  offset[n]++;
248  if(col_P[i] != n) offset[col_P[i]]++;
249  }
250  }
251  }
252 
253  // Divide the result by two
254  for(int i = 0; i < no_elem; i++) sym_val_P[i] /= 2.0;
255 
256  // Return symmetrized matrices
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;
260 
261  // Free up some memery
262  free(offset); offset = NULL;
263  free(row_counts); row_counts = NULL;
264  }
265 
266 private:
267 
268  void computeGradient(ScalarType* /*P*/, int* inp_row_P, int* inp_col_P, ScalarType* inp_val_P, ScalarType* Y, int N, int D, ScalarType* dC, ScalarType theta)
269  {
270  // Construct quadtree on current map
271  QuadTree* tree = new QuadTree(Y, N);
272 
273  // Compute all terms required for t-SNE gradient
274  ScalarType sum_Q = .0;
275  ScalarType* pos_f = (ScalarType*) calloc(N * D, sizeof(ScalarType));
276  ScalarType* neg_f = (ScalarType*) calloc(N * D, sizeof(ScalarType));
277  if(pos_f == NULL || neg_f == NULL) { printf("Memory allocation failed!\n"); exit(1); }
278  tree->computeEdgeForces(inp_row_P, inp_col_P, inp_val_P, N, pos_f);
279  for(int n = 0; n < N; n++) tree->computeNonEdgeForces(n, theta, neg_f + n * D, &sum_Q);
280 
281  // Compute final t-SNE gradient
282  for(int i = 0; i < N * D; i++) {
283  dC[i] = pos_f[i] - (neg_f[i] / sum_Q);
284  }
285  free(pos_f);
286  free(neg_f);
287  delete tree;
288  }
289 
290  void computeExactGradient(ScalarType* P, ScalarType* Y, int N, int D, ScalarType* dC)
291  {
292  // Make sure the current gradient contains zeros
293  for(int i = 0; i < N * D; i++) dC[i] = 0.0;
294 
295  // Compute the squared Euclidean distance matrix
296  ScalarType* DD = (ScalarType*) malloc(N * N * sizeof(ScalarType));
297  if(DD == NULL) { printf("Memory allocation failed!\n"); exit(1); }
298  computeSquaredEuclideanDistance(Y, N, D, DD);
299 
300  // Compute Q-matrix and normalization sum
301  ScalarType* Q = (ScalarType*) malloc(N * N * sizeof(ScalarType));
302  if(Q == NULL) { printf("Memory allocation failed!\n"); exit(1); }
303  ScalarType sum_Q = .0;
304  for(int n = 0; n < N; n++) {
305  for(int m = 0; m < N; m++) {
306  if(n != m) {
307  Q[n * N + m] = 1 / (1 + DD[n * N + m]);
308  sum_Q += Q[n * N + m];
309  }
310  }
311  }
312 
313  // Perform the computation of the gradient
314  for(int n = 0; n < N; n++) {
315  for(int m = 0; m < N; m++) {
316  if(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;
320  }
321  }
322  }
323  }
324 
325  // Free memory
326  free(DD); DD = NULL;
327  free(Q); Q = NULL;
328  }
329 
331  {
332  // Compute the squared Euclidean distance matrix
333  ScalarType* DD = (ScalarType*) malloc(N * N * sizeof(ScalarType));
334  ScalarType* Q = (ScalarType*) malloc(N * N * sizeof(ScalarType));
335  if(DD == NULL || Q == NULL) { printf("Memory allocation failed!\n"); exit(1); }
336  computeSquaredEuclideanDistance(Y, N, 2, DD);
337 
338  // Compute Q-matrix and normalization sum
339  ScalarType sum_Q = DBL_MIN;
340  for(int n = 0; n < N; n++) {
341  for(int m = 0; m < N; m++) {
342  if(n != m) {
343  Q[n * N + m] = 1 / (1 + DD[n * N + m]);
344  sum_Q += Q[n * N + m];
345  }
346  else Q[n * N + m] = DBL_MIN;
347  }
348  }
349  for(int i = 0; i < N * N; i++) Q[i] /= sum_Q;
350 
351  // Sum t-SNE error
352  ScalarType C = .0;
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));
356  }
357  }
358 
359  // Clean up memory
360  free(DD);
361  free(Q);
362  return C;
363  }
364 
365  ScalarType evaluateError(int* row_P, int* col_P, ScalarType* val_P, ScalarType* Y, int N, ScalarType theta)
366  {
367  // Get estimate of normalization term
368  const int QT_NO_DIMS = 2;
369  QuadTree* tree = new QuadTree(Y, N);
370  ScalarType buff[QT_NO_DIMS] = {.0, .0};
371  ScalarType sum_Q = .0;
372  for(int n = 0; n < N; n++) tree->computeNonEdgeForces(n, theta, buff, &sum_Q);
373 
374  // Loop over all edges to compute t-SNE error
375  int ind1, ind2;
376  ScalarType C = .0, Q;
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++) {
380  Q = .0;
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));
387  }
388  }
389  return C;
390  }
391 
392  void zeroMean(ScalarType* X, int N, int D)
393  {
394  // Compute data mean
395  ScalarType* mean = (ScalarType*) calloc(D, sizeof(ScalarType));
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];
400  }
401  }
402  for(int d = 0; d < D; d++) {
403  mean[d] /= (ScalarType) N;
404  }
405 
406  // Subtract data mean
407  for(int n = 0; n < N; n++) {
408  for(int d = 0; d < D; d++) {
409  X[n * D + d] -= mean[d];
410  }
411  }
412  free(mean); mean = NULL;
413  }
414 
415  void computeGaussianPerplexity(ScalarType* X, int N, int D, ScalarType* P, ScalarType perplexity)
416  {
417  // Compute the squared Euclidean distance matrix
418  ScalarType* DD = (ScalarType*) malloc(N * N * sizeof(ScalarType));
419  if(DD == NULL) { printf("Memory allocation failed!\n"); exit(1); }
420  computeSquaredEuclideanDistance(X, N, D, DD);
421 
422  // Compute the Gaussian kernel row by row
423  for(int n = 0; n < N; n++) {
424 
425  // Initialize some variables
426  bool found = false;
427  ScalarType beta = 1.0;
428  ScalarType min_beta = -DBL_MAX;
429  ScalarType max_beta = DBL_MAX;
430  ScalarType tol = 1e-5;
431  ScalarType sum_P;
432 
433  // Iterate until we found a good perplexity
434  int iter = 0;
435  while(!found && iter < 200) {
436 
437  // Compute Gaussian kernel row
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;
440 
441  // Compute entropy of current row
442  sum_P = DBL_MIN;
443  for(int m = 0; m < N; m++) sum_P += P[n * N + m];
444  ScalarType H = 0.0;
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);
447 
448  // Evaluate whether the entropy is within the tolerance level
449  ScalarType Hdiff = H - log(perplexity);
450  if(Hdiff < tol && -Hdiff < tol) {
451  found = true;
452  }
453  else {
454  if(Hdiff > 0) {
455  min_beta = beta;
456  if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
457  beta *= 2.0;
458  else
459  beta = (beta + max_beta) / 2.0;
460  }
461  else {
462  max_beta = beta;
463  if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
464  beta /= 2.0;
465  else
466  beta = (beta + min_beta) / 2.0;
467  }
468  }
469 
470  // Update iteration counter
471  iter++;
472  }
473 
474  // Row normalize P
475  for(int m = 0; m < N; m++) P[n * N + m] /= sum_P;
476  }
477 
478  // Clean up memory
479  free(DD); DD = NULL;
480  }
481 
482  void computeGaussianPerplexity(ScalarType* X, int N, int D, int** _row_P, int** _col_P, ScalarType** _val_P, ScalarType perplexity, int K)
483  {
484  if(perplexity > K) printf("Perplexity should be lower than K!\n");
485 
486  // Allocate the memory we need
487  *_row_P = (int*) malloc((N + 1) * sizeof(int));
488  *_col_P = (int*) calloc(N * K, sizeof(int));
489  *_val_P = (ScalarType*) calloc(N * K, sizeof(ScalarType));
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;
493  ScalarType* val_P = *_val_P;
494  ScalarType* cur_P = (ScalarType*) malloc((N - 1) * sizeof(ScalarType));
495  if(cur_P == NULL) { printf("Memory allocation failed!\n"); exit(1); }
496  row_P[0] = 0;
497  for(int n = 0; n < N; n++) row_P[n + 1] = row_P[n] + K;
498 
499  // Build ball tree on data set
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);
503  tree->create(obj_X);
504 
505  // Loop over all points to find nearest neighbors
506  //printf("Building tree...\n");
507  std::vector<DataPoint> indices;
508  std::vector<ScalarType> distances;
509  for(int n = 0; n < N; n++) {
510 
511  //if(n % 10000 == 0) printf(" - point %d of %d\n", n, N);
512 
513  // Find nearest neighbors
514  indices.clear();
515  distances.clear();
516  tree->search(obj_X[n], K + 1, &indices, &distances);
517 
518  // Initialize some variables for binary search
519  bool found = false;
520  ScalarType beta = 1.0;
521  ScalarType min_beta = -DBL_MAX;
522  ScalarType max_beta = DBL_MAX;
523  ScalarType tol = 1e-5;
524 
525  // Iterate until we found a good perplexity
526  int iter = 0; ScalarType sum_P;
527  while(!found && iter < 200) {
528 
529  // Compute Gaussian kernel row
530  for(int m = 0; m < K; m++) cur_P[m] = exp(-beta * distances[m + 1]);
531 
532  // Compute entropy of current row
533  sum_P = DBL_MIN;
534  for(int m = 0; m < K; m++) sum_P += cur_P[m];
535  ScalarType H = .0;
536  for(int m = 0; m < K; m++) H += beta * (distances[m + 1] * cur_P[m]);
537  H = (H / sum_P) + log(sum_P);
538 
539  // Evaluate whether the entropy is within the tolerance level
540  ScalarType Hdiff = H - log(perplexity);
541  if(Hdiff < tol && -Hdiff < tol) {
542  found = true;
543  }
544  else {
545  if(Hdiff > 0) {
546  min_beta = beta;
547  if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
548  beta *= 2.0;
549  else
550  beta = (beta + max_beta) / 2.0;
551  }
552  else {
553  max_beta = beta;
554  if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
555  beta /= 2.0;
556  else
557  beta = (beta + min_beta) / 2.0;
558  }
559  }
560 
561  // Update iteration counter
562  iter++;
563  }
564 
565  // Row-normalize current row of P and store in matrix
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];
570  }
571  }
572 
573  // Clean up memory
574  obj_X.clear();
575  free(cur_P);
576  delete tree;
577  }
578 
579  void computeGaussianPerplexity(ScalarType* X, int N, int D, int** _row_P, int** _col_P, ScalarType** _val_P, ScalarType perplexity, ScalarType threshold)
580  {
581  // Allocate some memory we need for computations
582  ScalarType* buff = (ScalarType*) malloc(D * sizeof(ScalarType));
583  ScalarType* DD = (ScalarType*) malloc(N * sizeof(ScalarType));
584  ScalarType* cur_P = (ScalarType*) malloc(N * sizeof(ScalarType));
585  if(buff == NULL || DD == NULL || cur_P == NULL) { printf("Memory allocation failed!\n"); exit(1); }
586 
587  // Compute the Gaussian kernel row by row (to find number of elements in sparse P)
588  int total_count = 0;
589  for(int n = 0; n < N; n++) {
590 
591  // Compute the squared Euclidean distance matrix
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];
595  DD[m] = .0;
596  for(int d = 0; d < D; d++) DD[m] += buff[d] * buff[d];
597  }
598 
599  // Initialize some variables
600  bool found = false;
601  ScalarType beta = 1.0;
602  ScalarType min_beta = -DBL_MAX;
603  ScalarType max_beta = DBL_MAX;
604  ScalarType tol = 1e-5;
605 
606  // Iterate until we found a good perplexity
607  int iter = 0; ScalarType sum_P;
608  while(!found && iter < 200) {
609 
610  // Compute Gaussian kernel row
611  for(int m = 0; m < N; m++) cur_P[m] = exp(-beta * DD[m]);
612  cur_P[n] = DBL_MIN;
613 
614  // Compute entropy of current row
615  sum_P = DBL_MIN;
616  for(int m = 0; m < N; m++) sum_P += cur_P[m];
617  ScalarType H = 0.0;
618  for(int m = 0; m < N; m++) H += beta * (DD[m] * cur_P[m]);
619  H = (H / sum_P) + log(sum_P);
620 
621  // Evaluate whether the entropy is within the tolerance level
622  ScalarType Hdiff = H - log(perplexity);
623  if(Hdiff < tol && -Hdiff < tol) {
624  found = true;
625  }
626  else {
627  if(Hdiff > 0) {
628  min_beta = beta;
629  if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
630  beta *= 2.0;
631  else
632  beta = (beta + max_beta) / 2.0;
633  }
634  else {
635  max_beta = beta;
636  if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
637  beta /= 2.0;
638  else
639  beta = (beta + min_beta) / 2.0;
640  }
641  }
642 
643  // Update iteration counter
644  iter++;
645  }
646 
647  // Row-normalize and threshold current row of P
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++;
651  }
652  }
653 
654  // Allocate the memory we need
655  *_row_P = (int*) malloc((N + 1) * sizeof(int));
656  *_col_P = (int*) malloc(total_count * sizeof(int));
657  *_val_P = (ScalarType*) malloc(total_count * sizeof(ScalarType));
658  int* row_P = *_row_P;
659  int* col_P = *_col_P;
660  ScalarType* val_P = *_val_P;
661  if(row_P == NULL || col_P == NULL || val_P == NULL) { printf("Memory allocation failed!\n"); exit(1); }
662  row_P[0] = 0;
663 
664  // Compute the Gaussian kernel row by row (this time, store the results)
665  int count = 0;
666  for(int n = 0; n < N; n++) {
667 
668  // Compute the squared Euclidean distance matrix
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];
672  DD[m] = .0;
673  for(int d = 0; d < D; d++) DD[m] += buff[d] * buff[d];
674  }
675 
676  // Initialize some variables
677  bool found = false;
678  ScalarType beta = 1.0;
679  ScalarType min_beta = -DBL_MAX;
680  ScalarType max_beta = DBL_MAX;
681  ScalarType tol = 1e-5;
682 
683  // Iterate until we found a good perplexity
684  int iter = 0; ScalarType sum_P;
685  while(!found && iter < 200) {
686 
687  // Compute Gaussian kernel row
688  for(int m = 0; m < N; m++) cur_P[m] = exp(-beta * DD[m]);
689  cur_P[n] = DBL_MIN;
690 
691  // Compute entropy of current row
692  sum_P = DBL_MIN;
693  for(int m = 0; m < N; m++) sum_P += cur_P[m];
694  ScalarType H = 0.0;
695  for(int m = 0; m < N; m++) H += beta * (DD[m] * cur_P[m]);
696  H = (H / sum_P) + log(sum_P);
697 
698  // Evaluate whether the entropy is within the tolerance level
699  ScalarType Hdiff = H - log(perplexity);
700  if(Hdiff < tol && -Hdiff < tol) {
701  found = true;
702  }
703  else {
704  if(Hdiff > 0) {
705  min_beta = beta;
706  if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
707  beta *= 2.0;
708  else
709  beta = (beta + max_beta) / 2.0;
710  }
711  else {
712  max_beta = beta;
713  if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
714  beta /= 2.0;
715  else
716  beta = (beta + min_beta) / 2.0;
717  }
718  }
719 
720  // Update iteration counter
721  iter++;
722  }
723 
724  // Row-normalize and threshold current row of P
725  for(int m = 0; m < N; m++) cur_P[m] /= sum_P;
726  for(int m = 0; m < N; m++) {
727  if(cur_P[m] > threshold / (ScalarType) N) {
728  col_P[count] = m;
729  val_P[count] = cur_P[m];
730  count++;
731  }
732  }
733  row_P[n + 1] = count;
734  }
735 
736  // Clean up memory
737  free(DD); DD = NULL;
738  free(buff); buff = NULL;
739  free(cur_P); cur_P = NULL;
740  }
741 
743  {
744  ScalarType* dataSums = (ScalarType*) calloc(N, sizeof(ScalarType));
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]);
749  }
750  }
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];
754  }
755  }
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;
759 
760  //cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans, N, N, D, -2.0, X, D, X, D, 1.0, DD, N);
761  free(dataSums); dataSums = NULL;
762  }
763 
764 };
765 
766 }
767 
768 #endif
Eigen::Matrix< tapkee::ScalarType, Eigen::Dynamic, Eigen::Dynamic > DenseMatrix
dense matrix type (non-overridable)
Definition: types.hpp:23
ScalarType evaluateError(ScalarType *P, ScalarType *Y, int N)
Definition: tsne.hpp:330
Namespace containing implementation of t-SNE algorithm.
Definition: quadtree.hpp:41
void computeExactGradient(ScalarType *P, ScalarType *Y, int N, int D, ScalarType *dC)
Definition: tsne.hpp:290
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) ...
Definition: types.hpp:15
void run(tapkee::DenseMatrix &X, int N, int D, ScalarType *Y, int no_dims, ScalarType perplexity, ScalarType theta)
Definition: tsne.hpp:60
void computeNonEdgeForces(int point_index, ScalarType theta, ScalarType neg_f[], ScalarType *sum_Q)
Definition: quadtree.hpp:317
std::string format(const std::string &fmt, const ValueWrapper a)
Definition: formatting.hpp:411
void computeGaussianPerplexity(ScalarType *X, int N, int D, int **_row_P, int **_col_P, ScalarType **_val_P, ScalarType perplexity, int K)
Definition: tsne.hpp:482
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)
Definition: tsne.hpp:268
void zeroMean(ScalarType *X, int N, int D)
Definition: tsne.hpp:392
void symmetrizeMatrix(int **_row_P, int **_col_P, ScalarType **_val_P, int N)
Definition: tsne.hpp:180
void computeSquaredEuclideanDistance(ScalarType *X, int N, int D, ScalarType *DD)
Definition: tsne.hpp:742
void create(const std::vector< T > &items)
void computeEdgeForces(int *row_P, int *col_P, ScalarType *val_P, int N, ScalarType *pos_f)
Definition: quadtree.hpp:350
void message_info(const std::string &msg)
Definition: logging.hpp:115
static LoggingSingleton & instance()
Definition: logging.hpp:102
ScalarType gaussian_random()
Definition: random.hpp:39
void computeGaussianPerplexity(ScalarType *X, int N, int D, int **_row_P, int **_col_P, ScalarType **_val_P, ScalarType perplexity, ScalarType threshold)
Definition: tsne.hpp:579
void computeGaussianPerplexity(ScalarType *X, int N, int D, ScalarType *P, ScalarType perplexity)
Definition: tsne.hpp:415
ScalarType evaluateError(int *row_P, int *col_P, ScalarType *val_P, ScalarType *Y, int N, ScalarType theta)
Definition: tsne.hpp:365
static ScalarType sign(ScalarType x)
Definition: tsne.hpp:55