Tapkee
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
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 static inline double sign(double x) { return (x == .0 ? .0 : (x < .0 ? -1.0 : 1.0)); }
54 
55 class TSNE
56 {
57 public:
58  void run(double* X, int N, int D, double* Y, int no_dims, double perplexity, double theta)
59  {
60  // Determine whether we are using an exact algorithm
61  bool exact = (theta == .0) ? true : false;
62  if (exact)
63  tapkee::LoggingSingleton::instance().message_info("Using exact t-SNE algorithm");
64  else
65  tapkee::LoggingSingleton::instance().message_info("Using Barnes-Hut-SNE algorithm");
66 
67  // Set learning parameters
68  float total_time = .0;
69  clock_t start, end;
70  int max_iter = 1000, stop_lying_iter = 250, mom_switch_iter = 250;
71  double momentum = .5, final_momentum = .8;
72  double eta = 200.0;
73 
74  // Allocate some memory
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;
81 
82  // Normalize input data (to prevent numerical problems)
83  double* P=NULL; int* row_P; int* col_P; double* val_P;
84  {
85  tapkee::tapkee_internal::timed_context context("Input similarities computation");
86  start = clock();
87  zeroMean(X, N, D);
88  double max_X = .0;
89  for(int i = 0; i < N * D; i++) {
90  if(X[i] > max_X) max_X = X[i];
91  }
92  for(int i = 0; i < N * D; i++) X[i] /= max_X;
93 
94  // Compute input similarities for exact t-SNE
95  if(exact) {
96 
97  // Compute similarities
98  P = (double*) malloc(N * N * sizeof(double));
99  if(P == NULL) { printf("Memory allocation failed!\n"); exit(1); }
100  computeGaussianPerplexity(X, N, D, P, perplexity);
101 
102  // Symmetrize input similarities
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];
107  }
108  }
109  double sum_P = .0;
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;
112  }
113 
114  // Compute input similarities for approximate t-SNE
115  else {
116 
117  // Compute asymmetric pairwise input similarities
118  computeGaussianPerplexity(X, N, D, &row_P, &col_P, &val_P, perplexity, (int) (3 * perplexity));
119 
120  // Symmetrize input similarities
121  symmetrizeMatrix(&row_P, &col_P, &val_P, N);
122  double sum_P = .0;
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;
125  }
126 
127  // Lie about the P-values
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; }
130 
131  // Initialize solution (randomly)
132  for(int i = 0; i < N * no_dims; i++) Y[i] = tapkee::gaussian_random() * .0001;
133  }
134 
135  {
136  tapkee::tapkee_internal::timed_context context("Main t-SNE loop");
137  for(int iter = 0; iter < max_iter; iter++) {
138 
139  // Compute (approximate) gradient
140  if(exact) computeExactGradient(P, Y, N, no_dims, dY);
141  else computeGradient(P, row_P, col_P, val_P, Y, N, no_dims, dY, theta);
142 
143  // Update gains
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;
146 
147  // Perform gradient update (with momentum and gains)
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];
150 
151  // Make solution zero-mean
152  zeroMean(Y, N, no_dims);
153 
154  // Stop lying about the P-values after a while, and switch momentum
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; }
158  }
159  if(iter == mom_switch_iter) momentum = final_momentum;
160 
161  // Print out progress
162  /*
163  if((iter > 0) && ((iter % 50 == 0) || (iter == max_iter - 1))) {
164  end = clock();
165  double C = .0;
166  if(exact) C = evaluateError(P, Y, N);
167  else C = evaluateError(row_P, col_P, val_P, Y, N, theta); // doing approximate computation here!
168  if(iter == 0)
169  {
170  //printf("Iteration %d: error is %f\n", iter + 1, C);
171  }
172  else
173  {
174  total_time += (float) (end - start) / CLOCKS_PER_SEC;
175  //printf("Iteration %d: error is %f (50 iterations in %4.2f seconds)\n", iter, C, (float) (end - start) / CLOCKS_PER_SEC);
176  }
177  start = clock();
178  }
179  */
180  }
181  end = clock(); total_time += (float) (end - start) / CLOCKS_PER_SEC;
182 
183  // Clean up memory
184  free(dY);
185  free(uY);
186  free(gains);
187  if(exact) free(P);
188  else {
189  free(row_P); row_P = NULL;
190  free(col_P); col_P = NULL;
191  free(val_P); val_P = NULL;
192  }
193  }
194  }
195 
196  void symmetrizeMatrix(int** _row_P, int** _col_P, double** _val_P, int N)
197  {
198  // Get sparse matrix
199  int* row_P = *_row_P;
200  int* col_P = *_col_P;
201  double* val_P = *_val_P;
202 
203  // Count number of elements and row counts of symmetric matrix
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++) {
208  // Check whether element (col_P[i], n) is present
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;
212  }
213  if(present) row_counts[n]++;
214  else {
215  row_counts[n]++;
216  row_counts[col_P[i]]++;
217  }
218  }
219  }
220  int no_elem = 0;
221  for(int n = 0; n < N; n++) no_elem += row_counts[n];
222 
223  // Allocate memory for symmetrized matrix
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); }
228 
229  // Construct new row indices for symmetric matrix
230  sym_row_P[0] = 0;
231  for(int n = 0; n < N; n++) sym_row_P[n + 1] = sym_row_P[n] + row_counts[n];
232 
233  // Fill the result matrix
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++) { // considering element(n, col_P[i])
238 
239  // Check whether element (col_P[i], n) is present
240  bool present = false;
241  for(int m = row_P[col_P[i]]; m < row_P[col_P[i] + 1]; m++) {
242  if(col_P[m] == n) {
243  present = true;
244  if(n <= col_P[i]) { // make sure we do not add elements twice
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];
249  }
250  }
251  }
252 
253  // If (col_P[i], n) is not present, there is no addition involved
254  if(!present) {
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];
259  }
260 
261  // Update offsets
262  if(!present || (present && n <= col_P[i])) {
263  offset[n]++;
264  if(col_P[i] != n) offset[col_P[i]]++;
265  }
266  }
267  }
268 
269  // Divide the result by two
270  for(int i = 0; i < no_elem; i++) sym_val_P[i] /= 2.0;
271 
272  // Return symmetrized matrices
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;
276 
277  // Free up some memery
278  free(offset); offset = NULL;
279  free(row_counts); row_counts = NULL;
280  }
281 
282 private:
283 
284  void computeGradient(double* /*P*/, int* inp_row_P, int* inp_col_P, double* inp_val_P, double* Y, int N, int D, double* dC, double theta)
285  {
286  // Construct quadtree on current map
287  QuadTree* tree = new QuadTree(Y, N);
288 
289  // Compute all terms required for t-SNE gradient
290  double sum_Q = .0;
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); }
294  tree->computeEdgeForces(inp_row_P, inp_col_P, inp_val_P, N, pos_f);
295  for(int n = 0; n < N; n++) tree->computeNonEdgeForces(n, theta, neg_f + n * D, &sum_Q);
296 
297  // Compute final t-SNE gradient
298  for(int i = 0; i < N * D; i++) {
299  dC[i] = pos_f[i] - (neg_f[i] / sum_Q);
300  }
301  free(pos_f);
302  free(neg_f);
303  delete tree;
304  }
305 
306  void computeExactGradient(double* P, double* Y, int N, int D, double* dC)
307  {
308  // Make sure the current gradient contains zeros
309  for(int i = 0; i < N * D; i++) dC[i] = 0.0;
310 
311  // Compute the squared Euclidean distance matrix
312  double* DD = (double*) malloc(N * N * sizeof(double));
313  if(DD == NULL) { printf("Memory allocation failed!\n"); exit(1); }
314  computeSquaredEuclideanDistance(Y, N, D, DD);
315 
316  // Compute Q-matrix and normalization sum
317  double* Q = (double*) malloc(N * N * sizeof(double));
318  if(Q == NULL) { printf("Memory allocation failed!\n"); exit(1); }
319  double sum_Q = .0;
320  for(int n = 0; n < N; n++) {
321  for(int m = 0; m < N; m++) {
322  if(n != m) {
323  Q[n * N + m] = 1 / (1 + DD[n * N + m]);
324  sum_Q += Q[n * N + m];
325  }
326  }
327  }
328 
329  // Perform the computation of the gradient
330  for(int n = 0; n < N; n++) {
331  for(int m = 0; m < N; m++) {
332  if(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;
336  }
337  }
338  }
339  }
340 
341  // Free memory
342  free(DD); DD = NULL;
343  free(Q); Q = NULL;
344  }
345 
346  double evaluateError(double* P, double* Y, int N)
347  {
348  // Compute the squared Euclidean distance matrix
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); }
352  computeSquaredEuclideanDistance(Y, N, 2, DD);
353 
354  // Compute Q-matrix and normalization sum
355  double sum_Q = DBL_MIN;
356  for(int n = 0; n < N; n++) {
357  for(int m = 0; m < N; m++) {
358  if(n != m) {
359  Q[n * N + m] = 1 / (1 + DD[n * N + m]);
360  sum_Q += Q[n * N + m];
361  }
362  else Q[n * N + m] = DBL_MIN;
363  }
364  }
365  for(int i = 0; i < N * N; i++) Q[i] /= sum_Q;
366 
367  // Sum t-SNE error
368  double C = .0;
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));
372  }
373  }
374 
375  // Clean up memory
376  free(DD);
377  free(Q);
378  return C;
379  }
380 
381  double evaluateError(int* row_P, int* col_P, double* val_P, double* Y, int N, double theta)
382  {
383  // Get estimate of normalization term
384  const int QT_NO_DIMS = 2;
385  QuadTree* tree = new QuadTree(Y, N);
386  double buff[QT_NO_DIMS] = {.0, .0};
387  double sum_Q = .0;
388  for(int n = 0; n < N; n++) tree->computeNonEdgeForces(n, theta, buff, &sum_Q);
389 
390  // Loop over all edges to compute t-SNE error
391  int ind1, ind2;
392  double C = .0, Q;
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++) {
396  Q = .0;
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));
403  }
404  }
405  return C;
406  }
407 
408  void zeroMean(double* X, int N, int D)
409  {
410  // Compute data mean
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];
416  }
417  }
418  for(int d = 0; d < D; d++) {
419  mean[d] /= (double) N;
420  }
421 
422  // Subtract data mean
423  for(int n = 0; n < N; n++) {
424  for(int d = 0; d < D; d++) {
425  X[n * D + d] -= mean[d];
426  }
427  }
428  free(mean); mean = NULL;
429  }
430 
431  void computeGaussianPerplexity(double* X, int N, int D, double* P, double perplexity)
432  {
433  // Compute the squared Euclidean distance matrix
434  double* DD = (double*) malloc(N * N * sizeof(double));
435  if(DD == NULL) { printf("Memory allocation failed!\n"); exit(1); }
436  computeSquaredEuclideanDistance(X, N, D, DD);
437 
438  // Compute the Gaussian kernel row by row
439  for(int n = 0; n < N; n++) {
440 
441  // Initialize some variables
442  bool found = false;
443  double beta = 1.0;
444  double min_beta = -DBL_MAX;
445  double max_beta = DBL_MAX;
446  double tol = 1e-5;
447  double sum_P;
448 
449  // Iterate until we found a good perplexity
450  int iter = 0;
451  while(!found && iter < 200) {
452 
453  // Compute Gaussian kernel row
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;
456 
457  // Compute entropy of current row
458  sum_P = DBL_MIN;
459  for(int m = 0; m < N; m++) sum_P += P[n * N + m];
460  double H = 0.0;
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);
463 
464  // Evaluate whether the entropy is within the tolerance level
465  double Hdiff = H - log(perplexity);
466  if(Hdiff < tol && -Hdiff < tol) {
467  found = true;
468  }
469  else {
470  if(Hdiff > 0) {
471  min_beta = beta;
472  if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
473  beta *= 2.0;
474  else
475  beta = (beta + max_beta) / 2.0;
476  }
477  else {
478  max_beta = beta;
479  if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
480  beta /= 2.0;
481  else
482  beta = (beta + min_beta) / 2.0;
483  }
484  }
485 
486  // Update iteration counter
487  iter++;
488  }
489 
490  // Row normalize P
491  for(int m = 0; m < N; m++) P[n * N + m] /= sum_P;
492  }
493 
494  // Clean up memory
495  free(DD); DD = NULL;
496  }
497 
498  void computeGaussianPerplexity(double* X, int N, int D, int** _row_P, int** _col_P, double** _val_P, double perplexity, int K)
499  {
500  if(perplexity > K) printf("Perplexity should be lower than K!\n");
501 
502  // Allocate the memory we need
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); }
512  row_P[0] = 0;
513  for(int n = 0; n < N; n++) row_P[n + 1] = row_P[n] + K;
514 
515  // Build ball tree on data set
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);
519  tree->create(obj_X);
520 
521  // Loop over all points to find nearest neighbors
522  //printf("Building tree...\n");
523  std::vector<DataPoint> indices;
524  std::vector<double> distances;
525  for(int n = 0; n < N; n++) {
526 
527  //if(n % 10000 == 0) printf(" - point %d of %d\n", n, N);
528 
529  // Find nearest neighbors
530  indices.clear();
531  distances.clear();
532  tree->search(obj_X[n], K + 1, &indices, &distances);
533 
534  // Initialize some variables for binary search
535  bool found = false;
536  double beta = 1.0;
537  double min_beta = -DBL_MAX;
538  double max_beta = DBL_MAX;
539  double tol = 1e-5;
540 
541  // Iterate until we found a good perplexity
542  int iter = 0; double sum_P;
543  while(!found && iter < 200) {
544 
545  // Compute Gaussian kernel row
546  for(int m = 0; m < K; m++) cur_P[m] = exp(-beta * distances[m + 1]);
547 
548  // Compute entropy of current row
549  sum_P = DBL_MIN;
550  for(int m = 0; m < K; m++) sum_P += cur_P[m];
551  double H = .0;
552  for(int m = 0; m < K; m++) H += beta * (distances[m + 1] * cur_P[m]);
553  H = (H / sum_P) + log(sum_P);
554 
555  // Evaluate whether the entropy is within the tolerance level
556  double Hdiff = H - log(perplexity);
557  if(Hdiff < tol && -Hdiff < tol) {
558  found = true;
559  }
560  else {
561  if(Hdiff > 0) {
562  min_beta = beta;
563  if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
564  beta *= 2.0;
565  else
566  beta = (beta + max_beta) / 2.0;
567  }
568  else {
569  max_beta = beta;
570  if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
571  beta /= 2.0;
572  else
573  beta = (beta + min_beta) / 2.0;
574  }
575  }
576 
577  // Update iteration counter
578  iter++;
579  }
580 
581  // Row-normalize current row of P and store in matrix
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];
586  }
587  }
588 
589  // Clean up memory
590  obj_X.clear();
591  free(cur_P);
592  delete tree;
593  }
594 
595  void computeGaussianPerplexity(double* X, int N, int D, int** _row_P, int** _col_P, double** _val_P, double perplexity, double threshold)
596  {
597  // Allocate some memory we need for computations
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); }
602 
603  // Compute the Gaussian kernel row by row (to find number of elements in sparse P)
604  int total_count = 0;
605  for(int n = 0; n < N; n++) {
606 
607  // Compute the squared Euclidean distance matrix
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];
611  DD[m] = .0;
612  for(int d = 0; d < D; d++) DD[m] += buff[d] * buff[d];
613  }
614 
615  // Initialize some variables
616  bool found = false;
617  double beta = 1.0;
618  double min_beta = -DBL_MAX;
619  double max_beta = DBL_MAX;
620  double tol = 1e-5;
621 
622  // Iterate until we found a good perplexity
623  int iter = 0; double sum_P;
624  while(!found && iter < 200) {
625 
626  // Compute Gaussian kernel row
627  for(int m = 0; m < N; m++) cur_P[m] = exp(-beta * DD[m]);
628  cur_P[n] = DBL_MIN;
629 
630  // Compute entropy of current row
631  sum_P = DBL_MIN;
632  for(int m = 0; m < N; m++) sum_P += cur_P[m];
633  double H = 0.0;
634  for(int m = 0; m < N; m++) H += beta * (DD[m] * cur_P[m]);
635  H = (H / sum_P) + log(sum_P);
636 
637  // Evaluate whether the entropy is within the tolerance level
638  double Hdiff = H - log(perplexity);
639  if(Hdiff < tol && -Hdiff < tol) {
640  found = true;
641  }
642  else {
643  if(Hdiff > 0) {
644  min_beta = beta;
645  if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
646  beta *= 2.0;
647  else
648  beta = (beta + max_beta) / 2.0;
649  }
650  else {
651  max_beta = beta;
652  if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
653  beta /= 2.0;
654  else
655  beta = (beta + min_beta) / 2.0;
656  }
657  }
658 
659  // Update iteration counter
660  iter++;
661  }
662 
663  // Row-normalize and threshold current row of P
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++;
667  }
668  }
669 
670  // Allocate the memory we need
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); }
678  row_P[0] = 0;
679 
680  // Compute the Gaussian kernel row by row (this time, store the results)
681  int count = 0;
682  for(int n = 0; n < N; n++) {
683 
684  // Compute the squared Euclidean distance matrix
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];
688  DD[m] = .0;
689  for(int d = 0; d < D; d++) DD[m] += buff[d] * buff[d];
690  }
691 
692  // Initialize some variables
693  bool found = false;
694  double beta = 1.0;
695  double min_beta = -DBL_MAX;
696  double max_beta = DBL_MAX;
697  double tol = 1e-5;
698 
699  // Iterate until we found a good perplexity
700  int iter = 0; double sum_P;
701  while(!found && iter < 200) {
702 
703  // Compute Gaussian kernel row
704  for(int m = 0; m < N; m++) cur_P[m] = exp(-beta * DD[m]);
705  cur_P[n] = DBL_MIN;
706 
707  // Compute entropy of current row
708  sum_P = DBL_MIN;
709  for(int m = 0; m < N; m++) sum_P += cur_P[m];
710  double H = 0.0;
711  for(int m = 0; m < N; m++) H += beta * (DD[m] * cur_P[m]);
712  H = (H / sum_P) + log(sum_P);
713 
714  // Evaluate whether the entropy is within the tolerance level
715  double Hdiff = H - log(perplexity);
716  if(Hdiff < tol && -Hdiff < tol) {
717  found = true;
718  }
719  else {
720  if(Hdiff > 0) {
721  min_beta = beta;
722  if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
723  beta *= 2.0;
724  else
725  beta = (beta + max_beta) / 2.0;
726  }
727  else {
728  max_beta = beta;
729  if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
730  beta /= 2.0;
731  else
732  beta = (beta + min_beta) / 2.0;
733  }
734  }
735 
736  // Update iteration counter
737  iter++;
738  }
739 
740  // Row-normalize and threshold current row of P
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) {
744  col_P[count] = m;
745  val_P[count] = cur_P[m];
746  count++;
747  }
748  }
749  row_P[n + 1] = count;
750  }
751 
752  // Clean up memory
753  free(DD); DD = NULL;
754  free(buff); buff = NULL;
755  free(cur_P); cur_P = NULL;
756  }
757 
758  void computeSquaredEuclideanDistance(double* X, int N, int D, double* DD)
759  {
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]);
765  }
766  }
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];
770  }
771  }
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;
775 
776  //cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans, N, N, D, -2.0, X, D, X, D, 1.0, DD, N);
777  free(dataSums); dataSums = NULL;
778  }
779 
780 };
781 
782 }
783 
784 #endif
void run(double *X, int N, int D, double *Y, int no_dims, double perplexity, double theta)
Definition: tsne.hpp:58
void symmetrizeMatrix(int **_row_P, int **_col_P, double **_val_P, int N)
Definition: tsne.hpp:196
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)
Definition: tsne.hpp:498
double evaluateError(double *P, double *Y, int N)
Definition: tsne.hpp:346
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)
Definition: tsne.hpp:284
double evaluateError(int *row_P, int *col_P, double *val_P, double *Y, int N, double theta)
Definition: tsne.hpp:381
void computeEdgeForces(int *row_P, int *col_P, double *val_P, int N, double *pos_f)
Definition: quadtree.hpp:348
void computeGaussianPerplexity(double *X, int N, int D, int **_row_P, int **_col_P, double **_val_P, double perplexity, double threshold)
Definition: tsne.hpp:595
static double sign(double x)
Definition: tsne.hpp:53
void create(const std::vector< T > &items)
void computeGaussianPerplexity(double *X, int N, int D, double *P, double perplexity)
Definition: tsne.hpp:431
void message_info(const std::string &msg)
Definition: logging.hpp:113
static LoggingSingleton & instance()
Definition: logging.hpp:100
ScalarType gaussian_random()
Definition: random.hpp:39
void computeExactGradient(double *P, double *Y, int N, int D, double *dC)
Definition: tsne.hpp:306
void computeSquaredEuclideanDistance(double *X, int N, int D, double *DD)
Definition: tsne.hpp:758
void zeroMean(double *X, int N, int D)
Definition: tsne.hpp:408
void computeNonEdgeForces(int point_index, double theta, double neg_f[], double *sum_Q)
Definition: quadtree.hpp:315