Tapkee
tsne.hpp
Go to the documentation of this file.
00001 
00032 #ifndef TSNE_H
00033 #define TSNE_H
00034 
00035 /* Tapkee includes */
00036 #include <tapkee/utils/logging.hpp>
00037 #include <tapkee/utils/time.hpp>
00038 #include <tapkee/external/barnes_hut_sne/quadtree.hpp>
00039 #include <tapkee/external/barnes_hut_sne/vptree.hpp>
00040 /* End of Tapkee includes */
00041 
00042 #include <math.h>
00043 #include <float.h>
00044 #include <stdlib.h>
00045 #include <stdio.h>
00046 #include <cstring>
00047 #include <time.h>
00048 
00050 namespace tsne
00051 {
00052 
00053 static inline double sign(double x) { return (x == .0 ? .0 : (x < .0 ? -1.0 : 1.0)); }
00054 
00055 class TSNE
00056 {    
00057 public:
00058     void run(double* X, int N, int D, double* Y, int no_dims, double perplexity, double theta)
00059     {
00060         // Determine whether we are using an exact algorithm
00061         bool exact = (theta == .0) ? true : false;
00062         if (exact) 
00063             tapkee::LoggingSingleton::instance().message_info("Using exact t-SNE algorithm");
00064         else
00065             tapkee::LoggingSingleton::instance().message_info("Using Barnes-Hut-SNE algorithm");
00066         
00067         // Set learning parameters
00068         float total_time = .0;
00069         clock_t start, end;
00070         int max_iter = 1000, stop_lying_iter = 250, mom_switch_iter = 250;
00071         double momentum = .5, final_momentum = .8;
00072         double eta = 200.0;
00073         
00074         // Allocate some memory
00075         double* dY    = (double*) malloc(N * no_dims * sizeof(double));
00076         double* uY    = (double*) malloc(N * no_dims * sizeof(double));
00077         double* gains = (double*) malloc(N * no_dims * sizeof(double));
00078         if(dY == NULL || uY == NULL || gains == NULL) { printf("Memory allocation failed!\n"); exit(1); }
00079         for(int i = 0; i < N * no_dims; i++)    uY[i] =  .0;
00080         for(int i = 0; i < N * no_dims; i++) gains[i] = 1.0;
00081         
00082         // Normalize input data (to prevent numerical problems)
00083         double* P=NULL; int* row_P; int* col_P; double* val_P;
00084         {
00085             tapkee::tapkee_internal::timed_context context("Input similarities computation");
00086             start = clock();
00087             zeroMean(X, N, D);
00088             double max_X = .0;
00089             for(int i = 0; i < N * D; i++) {
00090                 if(X[i] > max_X) max_X = X[i];
00091             }
00092             for(int i = 0; i < N * D; i++) X[i] /= max_X;
00093             
00094             // Compute input similarities for exact t-SNE
00095             if(exact) {
00096                 
00097                 // Compute similarities
00098                 P = (double*) malloc(N * N * sizeof(double));
00099                 if(P == NULL) { printf("Memory allocation failed!\n"); exit(1); }
00100                 computeGaussianPerplexity(X, N, D, P, perplexity);
00101             
00102                 // Symmetrize input similarities
00103                 for(int n = 0; n < N; n++) {
00104                     for(int m = n + 1; m < N; m++) {
00105                         P[n * N + m] += P[m * N + n];
00106                         P[m * N + n]  = P[n * N + m];
00107                     }
00108                 }
00109                 double sum_P = .0;
00110                 for(int i = 0; i < N * N; i++) sum_P += P[i];
00111                 for(int i = 0; i < N * N; i++) P[i] /= sum_P;
00112             }
00113             
00114             // Compute input similarities for approximate t-SNE
00115             else {
00116             
00117                 // Compute asymmetric pairwise input similarities
00118                 computeGaussianPerplexity(X, N, D, &row_P, &col_P, &val_P, perplexity, (int) (3 * perplexity));
00119                 
00120                 // Symmetrize input similarities
00121                 symmetrizeMatrix(&row_P, &col_P, &val_P, N);
00122                 double sum_P = .0;
00123                 for(int i = 0; i < row_P[N]; i++) sum_P += val_P[i];
00124                 for(int i = 0; i < row_P[N]; i++) val_P[i] /= sum_P;
00125             }
00126             
00127             // Lie about the P-values
00128             if(exact) { for(int i = 0; i < N * N; i++)        P[i] *= 12.0; }
00129             else {      for(int i = 0; i < row_P[N]; i++) val_P[i] *= 12.0; }
00130 
00131             // Initialize solution (randomly)
00132             for(int i = 0; i < N * no_dims; i++) Y[i] = tapkee::gaussian_random() * .0001;
00133         }
00134 
00135         {
00136             tapkee::tapkee_internal::timed_context context("Main t-SNE loop");
00137             for(int iter = 0; iter < max_iter; iter++) {
00138                 
00139                 // Compute (approximate) gradient
00140                 if(exact) computeExactGradient(P, Y, N, no_dims, dY);
00141                 else computeGradient(P, row_P, col_P, val_P, Y, N, no_dims, dY, theta);
00142                 
00143                 // Update gains
00144                 for(int i = 0; i < N * no_dims; i++) gains[i] = (sign(dY[i]) != sign(uY[i])) ? (gains[i] + .2) : (gains[i] * .8);
00145                 for(int i = 0; i < N * no_dims; i++) if(gains[i] < .01) gains[i] = .01;
00146                     
00147                 // Perform gradient update (with momentum and gains)
00148                 for(int i = 0; i < N * no_dims; i++) uY[i] = momentum * uY[i] - eta * gains[i] * dY[i];
00149                 for(int i = 0; i < N * no_dims; i++)  Y[i] = Y[i] + uY[i];
00150                 
00151                 // Make solution zero-mean
00152                 zeroMean(Y, N, no_dims);
00153                 
00154                 // Stop lying about the P-values after a while, and switch momentum
00155                 if(iter == stop_lying_iter) {
00156                     if(exact) { for(int i = 0; i < N * N; i++)        P[i] /= 12.0; }
00157                     else      { for(int i = 0; i < row_P[N]; i++) val_P[i] /= 12.0; }
00158                 }
00159                 if(iter == mom_switch_iter) momentum = final_momentum;
00160                 
00161                 // Print out progress
00162                 /*
00163                 if((iter > 0) && ((iter % 50 == 0) || (iter == max_iter - 1))) {
00164                     end = clock();
00165                     double C = .0;
00166                     if(exact) C = evaluateError(P, Y, N);
00167                     else      C = evaluateError(row_P, col_P, val_P, Y, N, theta);  // doing approximate computation here!
00168                     if(iter == 0) 
00169                     {
00170                         //printf("Iteration %d: error is %f\n", iter + 1, C);
00171                     }
00172                     else 
00173                     {
00174                         total_time += (float) (end - start) / CLOCKS_PER_SEC;
00175                         //printf("Iteration %d: error is %f (50 iterations in %4.2f seconds)\n", iter, C, (float) (end - start) / CLOCKS_PER_SEC);
00176                     }
00177                     start = clock();
00178                 }
00179                 */
00180             }
00181             end = clock(); total_time += (float) (end - start) / CLOCKS_PER_SEC;
00182             
00183             // Clean up memory
00184             free(dY);
00185             free(uY);
00186             free(gains);
00187             if(exact) free(P);
00188             else {
00189                 free(row_P); row_P = NULL;
00190                 free(col_P); col_P = NULL;
00191                 free(val_P); val_P = NULL;
00192             }
00193         }
00194     }
00195 
00196     void symmetrizeMatrix(int** _row_P, int** _col_P, double** _val_P, int N)
00197     { 
00198         // Get sparse matrix
00199         int* row_P = *_row_P;
00200         int* col_P = *_col_P;
00201         double* val_P = *_val_P;
00202 
00203         // Count number of elements and row counts of symmetric matrix
00204         int* row_counts = (int*) calloc(N, sizeof(int));
00205         if(row_counts == NULL) { printf("Memory allocation failed!\n"); exit(1); }
00206         for(int n = 0; n < N; n++) {
00207             for(int i = row_P[n]; i < row_P[n + 1]; i++) {
00208                 // Check whether element (col_P[i], n) is present
00209                 bool present = false;
00210                 for(int m = row_P[col_P[i]]; m < row_P[col_P[i] + 1]; m++) {
00211                     if(col_P[m] == n) present = true;
00212                 }
00213                 if(present) row_counts[n]++;
00214                 else {
00215                     row_counts[n]++;
00216                     row_counts[col_P[i]]++;
00217                 }
00218             }
00219         }
00220         int no_elem = 0;
00221         for(int n = 0; n < N; n++) no_elem += row_counts[n];
00222         
00223         // Allocate memory for symmetrized matrix
00224         int*    sym_row_P = (int*)    malloc((N + 1) * sizeof(int));
00225         int*    sym_col_P = (int*)    malloc(no_elem * sizeof(int));
00226         double* sym_val_P = (double*) malloc(no_elem * sizeof(double));
00227         if(sym_row_P == NULL || sym_col_P == NULL || sym_val_P == NULL) { printf("Memory allocation failed!\n"); exit(1); }
00228         
00229         // Construct new row indices for symmetric matrix
00230         sym_row_P[0] = 0;
00231         for(int n = 0; n < N; n++) sym_row_P[n + 1] = sym_row_P[n] + row_counts[n];
00232         
00233         // Fill the result matrix
00234         int* offset = (int*) calloc(N, sizeof(int));
00235         if(offset == NULL) { printf("Memory allocation failed!\n"); exit(1); }
00236         for(int n = 0; n < N; n++) {
00237             for(int i = row_P[n]; i < row_P[n + 1]; i++) {                                  // considering element(n, col_P[i])
00238                 
00239                 // Check whether element (col_P[i], n) is present
00240                 bool present = false;
00241                 for(int m = row_P[col_P[i]]; m < row_P[col_P[i] + 1]; m++) {
00242                     if(col_P[m] == n) {
00243                         present = true;
00244                         if(n <= col_P[i]) {                                                 // make sure we do not add elements twice
00245                             sym_col_P[sym_row_P[n]        + offset[n]]        = col_P[i];
00246                             sym_col_P[sym_row_P[col_P[i]] + offset[col_P[i]]] = n;
00247                             sym_val_P[sym_row_P[n]        + offset[n]]        = val_P[i] + val_P[m];
00248                             sym_val_P[sym_row_P[col_P[i]] + offset[col_P[i]]] = val_P[i] + val_P[m];
00249                         }
00250                     }
00251                 }
00252                 
00253                 // If (col_P[i], n) is not present, there is no addition involved
00254                 if(!present) {
00255                     sym_col_P[sym_row_P[n]        + offset[n]]        = col_P[i];
00256                     sym_col_P[sym_row_P[col_P[i]] + offset[col_P[i]]] = n;
00257                     sym_val_P[sym_row_P[n]        + offset[n]]        = val_P[i];
00258                     sym_val_P[sym_row_P[col_P[i]] + offset[col_P[i]]] = val_P[i];
00259                 }
00260                 
00261                 // Update offsets
00262                 if(!present || (present && n <= col_P[i])) {
00263                     offset[n]++;
00264                     if(col_P[i] != n) offset[col_P[i]]++;               
00265                 }
00266             }
00267         }
00268         
00269         // Divide the result by two
00270         for(int i = 0; i < no_elem; i++) sym_val_P[i] /= 2.0;
00271         
00272         // Return symmetrized matrices
00273         free(*_row_P); *_row_P = sym_row_P;
00274         free(*_col_P); *_col_P = sym_col_P;
00275         free(*_val_P); *_val_P = sym_val_P;
00276         
00277         // Free up some memery
00278         free(offset); offset = NULL;
00279         free(row_counts); row_counts  = NULL;
00280     }
00281     
00282 private:
00283 
00284     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)
00285     {
00286         // Construct quadtree on current map
00287         QuadTree* tree = new QuadTree(Y, N);
00288         
00289         // Compute all terms required for t-SNE gradient
00290         double sum_Q = .0;
00291         double* pos_f = (double*) calloc(N * D, sizeof(double));
00292         double* neg_f = (double*) calloc(N * D, sizeof(double));
00293         if(pos_f == NULL || neg_f == NULL) { printf("Memory allocation failed!\n"); exit(1); }
00294         tree->computeEdgeForces(inp_row_P, inp_col_P, inp_val_P, N, pos_f);
00295         for(int n = 0; n < N; n++) tree->computeNonEdgeForces(n, theta, neg_f + n * D, &sum_Q);
00296         
00297         // Compute final t-SNE gradient
00298         for(int i = 0; i < N * D; i++) {
00299             dC[i] = pos_f[i] - (neg_f[i] / sum_Q);
00300         }
00301         free(pos_f);
00302         free(neg_f);
00303         delete tree;
00304     }
00305 
00306     void computeExactGradient(double* P, double* Y, int N, int D, double* dC)
00307     {
00308         // Make sure the current gradient contains zeros
00309         for(int i = 0; i < N * D; i++) dC[i] = 0.0;
00310         
00311         // Compute the squared Euclidean distance matrix
00312         double* DD = (double*) malloc(N * N * sizeof(double));
00313         if(DD == NULL) { printf("Memory allocation failed!\n"); exit(1); }
00314         computeSquaredEuclideanDistance(Y, N, D, DD);
00315         
00316         // Compute Q-matrix and normalization sum
00317         double* Q    = (double*) malloc(N * N * sizeof(double));
00318         if(Q == NULL) { printf("Memory allocation failed!\n"); exit(1); }
00319         double sum_Q = .0;
00320         for(int n = 0; n < N; n++) {
00321             for(int m = 0; m < N; m++) {
00322                 if(n != m) {
00323                     Q[n * N + m] = 1 / (1 + DD[n * N + m]);
00324                     sum_Q += Q[n * N + m];
00325                 }
00326             }
00327         }
00328         
00329         // Perform the computation of the gradient
00330         for(int n = 0; n < N; n++) {
00331             for(int m = 0; m < N; m++) {
00332                 if(n != m) {
00333                     double mult = (P[n * N + m] - (Q[n * N + m] / sum_Q)) * Q[n * N + m];
00334                     for(int d = 0; d < D; d++) {
00335                         dC[n * D + d] += (Y[n * D + d] - Y[m * D + d]) * mult;
00336                     }
00337                 }
00338             }
00339         }
00340         
00341         // Free memory
00342         free(DD); DD = NULL;
00343         free(Q);  Q  = NULL;
00344     }
00345 
00346     double evaluateError(double* P, double* Y, int N)
00347     { 
00348         // Compute the squared Euclidean distance matrix
00349         double* DD = (double*) malloc(N * N * sizeof(double));
00350         double* Q = (double*) malloc(N * N * sizeof(double));
00351         if(DD == NULL || Q == NULL) { printf("Memory allocation failed!\n"); exit(1); }
00352         computeSquaredEuclideanDistance(Y, N, 2, DD);
00353         
00354         // Compute Q-matrix and normalization sum
00355         double sum_Q = DBL_MIN;
00356         for(int n = 0; n < N; n++) {
00357             for(int m = 0; m < N; m++) {
00358                 if(n != m) {
00359                     Q[n * N + m] = 1 / (1 + DD[n * N + m]);
00360                     sum_Q += Q[n * N + m];
00361                 }
00362                 else Q[n * N + m] = DBL_MIN;
00363             }
00364         }
00365         for(int i = 0; i < N * N; i++) Q[i] /= sum_Q;
00366         
00367         // Sum t-SNE error
00368         double C = .0;
00369         for(int n = 0; n < N; n++) {
00370             for(int m = 0; m < N; m++) {
00371                 C += P[n * N + m] * log((P[n * N + m] + 1e-9) / (Q[n * N + m] + 1e-9));
00372             }
00373         }
00374         
00375         // Clean up memory
00376         free(DD);
00377         free(Q);
00378         return C;
00379     }
00380 
00381     double evaluateError(int* row_P, int* col_P, double* val_P, double* Y, int N, double theta)
00382     {
00383         // Get estimate of normalization term
00384         const int QT_NO_DIMS = 2;
00385         QuadTree* tree = new QuadTree(Y, N);
00386         double buff[QT_NO_DIMS] = {.0, .0};
00387         double sum_Q = .0;
00388         for(int n = 0; n < N; n++) tree->computeNonEdgeForces(n, theta, buff, &sum_Q);
00389         
00390         // Loop over all edges to compute t-SNE error
00391         int ind1, ind2;
00392         double C = .0, Q;
00393         for(int n = 0; n < N; n++) {
00394             ind1 = n * QT_NO_DIMS;
00395             for(int i = row_P[n]; i < row_P[n + 1]; i++) {
00396                 Q = .0;
00397                 ind2 = col_P[i] * QT_NO_DIMS;
00398                 for(int d = 0; d < QT_NO_DIMS; d++) buff[d]  = Y[ind1 + d];
00399                 for(int d = 0; d < QT_NO_DIMS; d++) buff[d] -= Y[ind2 + d];
00400                 for(int d = 0; d < QT_NO_DIMS; d++) Q += buff[d] * buff[d];
00401                 Q = (1.0 / (1.0 + Q)) / sum_Q;
00402                 C += val_P[i] * log((val_P[i] + FLT_MIN) / (Q + FLT_MIN));
00403             }
00404         }
00405         return C;
00406     }
00407 
00408     void zeroMean(double* X, int N, int D)
00409     {
00410         // Compute data mean
00411         double* mean = (double*) calloc(D, sizeof(double));
00412         if(mean == NULL) { printf("Memory allocation failed!\n"); exit(1); }
00413         for(int n = 0; n < N; n++) {
00414             for(int d = 0; d < D; d++) {
00415                 mean[d] += X[n * D + d];
00416             }
00417         }
00418         for(int d = 0; d < D; d++) {
00419             mean[d] /= (double) N;
00420         }
00421         
00422         // Subtract data mean
00423         for(int n = 0; n < N; n++) {
00424             for(int d = 0; d < D; d++) {
00425                 X[n * D + d] -= mean[d];
00426             }
00427         }
00428         free(mean); mean = NULL;
00429     }
00430 
00431     void computeGaussianPerplexity(double* X, int N, int D, double* P, double perplexity)
00432     {
00433         // Compute the squared Euclidean distance matrix
00434         double* DD = (double*) malloc(N * N * sizeof(double));
00435         if(DD == NULL) { printf("Memory allocation failed!\n"); exit(1); }
00436         computeSquaredEuclideanDistance(X, N, D, DD);
00437         
00438         // Compute the Gaussian kernel row by row
00439         for(int n = 0; n < N; n++) {
00440             
00441             // Initialize some variables
00442             bool found = false;
00443             double beta = 1.0;
00444             double min_beta = -DBL_MAX;
00445             double max_beta =  DBL_MAX;
00446             double tol = 1e-5;
00447             double sum_P;
00448             
00449             // Iterate until we found a good perplexity
00450             int iter = 0;
00451             while(!found && iter < 200) {
00452                 
00453                 // Compute Gaussian kernel row
00454                 for(int m = 0; m < N; m++) P[n * N + m] = exp(-beta * DD[n * N + m]);
00455                 P[n * N + n] = DBL_MIN;
00456                 
00457                 // Compute entropy of current row
00458                 sum_P = DBL_MIN;
00459                 for(int m = 0; m < N; m++) sum_P += P[n * N + m];
00460                 double H = 0.0;
00461                 for(int m = 0; m < N; m++) H += beta * (DD[n * N + m] * P[n * N + m]);
00462                 H = (H / sum_P) + log(sum_P);
00463                 
00464                 // Evaluate whether the entropy is within the tolerance level
00465                 double Hdiff = H - log(perplexity);
00466                 if(Hdiff < tol && -Hdiff < tol) {
00467                     found = true;
00468                 }
00469                 else {
00470                     if(Hdiff > 0) {
00471                         min_beta = beta;
00472                         if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
00473                             beta *= 2.0;
00474                         else
00475                             beta = (beta + max_beta) / 2.0;
00476                     }
00477                     else {
00478                         max_beta = beta;
00479                         if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
00480                             beta /= 2.0;
00481                         else
00482                             beta = (beta + min_beta) / 2.0;
00483                     }
00484                 }
00485                 
00486                 // Update iteration counter
00487                 iter++;
00488             }
00489             
00490             // Row normalize P
00491             for(int m = 0; m < N; m++) P[n * N + m] /= sum_P;
00492         }
00493         
00494         // Clean up memory
00495         free(DD); DD = NULL;
00496     }
00497 
00498     void computeGaussianPerplexity(double* X, int N, int D, int** _row_P, int** _col_P, double** _val_P, double perplexity, int K)
00499     { 
00500         if(perplexity > K) printf("Perplexity should be lower than K!\n");
00501         
00502         // Allocate the memory we need
00503         *_row_P = (int*)    malloc((N + 1) * sizeof(int));
00504         *_col_P = (int*)    calloc(N * K, sizeof(int));
00505         *_val_P = (double*) calloc(N * K, sizeof(double));
00506         if(*_row_P == NULL || *_col_P == NULL || *_val_P == NULL) { printf("Memory allocation failed!\n"); exit(1); }
00507         int* row_P = *_row_P;
00508         int* col_P = *_col_P;
00509         double* val_P = *_val_P;
00510         double* cur_P = (double*) malloc((N - 1) * sizeof(double));
00511         if(cur_P == NULL) { printf("Memory allocation failed!\n"); exit(1); }
00512         row_P[0] = 0;
00513         for(int n = 0; n < N; n++) row_P[n + 1] = row_P[n] + K;    
00514         
00515         // Build ball tree on data set
00516         VpTree<DataPoint, euclidean_distance>* tree = new VpTree<DataPoint, euclidean_distance>();
00517         std::vector<DataPoint> obj_X(N, DataPoint(D, -1, X));
00518         for(int n = 0; n < N; n++) obj_X[n] = DataPoint(D, n, X + n * D);
00519         tree->create(obj_X);
00520         
00521         // Loop over all points to find nearest neighbors
00522         //printf("Building tree...\n");
00523         std::vector<DataPoint> indices;
00524         std::vector<double> distances;
00525         for(int n = 0; n < N; n++) {
00526             
00527             //if(n % 10000 == 0) printf(" - point %d of %d\n", n, N);
00528             
00529             // Find nearest neighbors
00530             indices.clear();
00531             distances.clear();
00532             tree->search(obj_X[n], K + 1, &indices, &distances);
00533             
00534             // Initialize some variables for binary search
00535             bool found = false;
00536             double beta = 1.0;
00537             double min_beta = -DBL_MAX;
00538             double max_beta =  DBL_MAX;
00539             double tol = 1e-5;
00540             
00541             // Iterate until we found a good perplexity
00542             int iter = 0; double sum_P;
00543             while(!found && iter < 200) {
00544                 
00545                 // Compute Gaussian kernel row
00546                 for(int m = 0; m < K; m++) cur_P[m] = exp(-beta * distances[m + 1]);
00547                 
00548                 // Compute entropy of current row
00549                 sum_P = DBL_MIN;
00550                 for(int m = 0; m < K; m++) sum_P += cur_P[m];
00551                 double H = .0;
00552                 for(int m = 0; m < K; m++) H += beta * (distances[m + 1] * cur_P[m]);
00553                 H = (H / sum_P) + log(sum_P);
00554                 
00555                 // Evaluate whether the entropy is within the tolerance level
00556                 double Hdiff = H - log(perplexity);
00557                 if(Hdiff < tol && -Hdiff < tol) {
00558                     found = true;
00559                 }
00560                 else {
00561                     if(Hdiff > 0) {
00562                         min_beta = beta;
00563                         if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
00564                             beta *= 2.0;
00565                         else
00566                             beta = (beta + max_beta) / 2.0;
00567                     }
00568                     else {
00569                         max_beta = beta;
00570                         if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
00571                             beta /= 2.0;
00572                         else
00573                             beta = (beta + min_beta) / 2.0;
00574                     }
00575                 }
00576                 
00577                 // Update iteration counter
00578                 iter++;
00579             }
00580             
00581             // Row-normalize current row of P and store in matrix
00582             for(int m = 0; m < K; m++) cur_P[m] /= sum_P;
00583             for(int m = 0; m < K; m++) {
00584                 col_P[row_P[n] + m] = indices[m + 1].index();
00585                 val_P[row_P[n] + m] = cur_P[m];
00586             }
00587         }
00588         
00589         // Clean up memory
00590         obj_X.clear();
00591         free(cur_P);
00592         delete tree;
00593     }
00594 
00595     void computeGaussianPerplexity(double* X, int N, int D, int** _row_P, int** _col_P, double** _val_P, double perplexity, double threshold)
00596     {
00597         // Allocate some memory we need for computations
00598         double* buff  = (double*) malloc(D * sizeof(double));
00599         double* DD    = (double*) malloc(N * sizeof(double));
00600         double* cur_P = (double*) malloc(N * sizeof(double));
00601         if(buff == NULL || DD == NULL || cur_P == NULL) { printf("Memory allocation failed!\n"); exit(1); }
00602 
00603         // Compute the Gaussian kernel row by row (to find number of elements in sparse P)
00604         int total_count = 0;
00605         for(int n = 0; n < N; n++) {
00606         
00607             // Compute the squared Euclidean distance matrix
00608             for(int m = 0; m < N; m++) {
00609                 for(int d = 0; d < D; d++) buff[d]  = X[n * D + d];
00610                 for(int d = 0; d < D; d++) buff[d] -= X[m * D + d];
00611                 DD[m] = .0;
00612                 for(int d = 0; d < D; d++) DD[m] += buff[d] * buff[d];
00613             }
00614            
00615             // Initialize some variables
00616             bool found = false;
00617             double beta = 1.0;
00618             double min_beta = -DBL_MAX;
00619             double max_beta =  DBL_MAX;
00620             double tol = 1e-5;
00621             
00622             // Iterate until we found a good perplexity
00623             int iter = 0; double sum_P;
00624             while(!found && iter < 200) {
00625                 
00626                 // Compute Gaussian kernel row
00627                 for(int m = 0; m < N; m++) cur_P[m] = exp(-beta * DD[m]);
00628                 cur_P[n] = DBL_MIN;
00629                 
00630                 // Compute entropy of current row
00631                 sum_P = DBL_MIN;
00632                 for(int m = 0; m < N; m++) sum_P += cur_P[m];
00633                 double H = 0.0;
00634                 for(int m = 0; m < N; m++) H += beta * (DD[m] * cur_P[m]);
00635                 H = (H / sum_P) + log(sum_P);
00636                 
00637                 // Evaluate whether the entropy is within the tolerance level
00638                 double Hdiff = H - log(perplexity);
00639                 if(Hdiff < tol && -Hdiff < tol) {
00640                     found = true;
00641                 }
00642                 else {
00643                     if(Hdiff > 0) {
00644                         min_beta = beta;
00645                         if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
00646                             beta *= 2.0;
00647                         else
00648                             beta = (beta + max_beta) / 2.0;
00649                     }
00650                     else {
00651                         max_beta = beta;
00652                         if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
00653                             beta /= 2.0;
00654                         else
00655                             beta = (beta + min_beta) / 2.0;
00656                     }
00657                 }
00658                 
00659                 // Update iteration counter
00660                 iter++;
00661             }
00662             
00663             // Row-normalize and threshold current row of P
00664             for(int m = 0; m < N; m++) cur_P[m] /= sum_P;
00665             for(int m = 0; m < N; m++) {
00666                 if(cur_P[m] > threshold / (double) N) total_count++;
00667             }
00668         }
00669         
00670         // Allocate the memory we need
00671         *_row_P = (int*)    malloc((N + 1)     * sizeof(int));
00672         *_col_P = (int*)    malloc(total_count * sizeof(int));
00673         *_val_P = (double*) malloc(total_count * sizeof(double));
00674         int* row_P = *_row_P;
00675         int* col_P = *_col_P;
00676         double* val_P = *_val_P;
00677         if(row_P == NULL || col_P == NULL || val_P == NULL) { printf("Memory allocation failed!\n"); exit(1); }
00678         row_P[0] = 0;
00679         
00680         // Compute the Gaussian kernel row by row (this time, store the results)
00681         int count = 0;
00682         for(int n = 0; n < N; n++) {
00683             
00684             // Compute the squared Euclidean distance matrix
00685             for(int m = 0; m < N; m++) {
00686                 for(int d = 0; d < D; d++) buff[d]  = X[n * D + d];
00687                 for(int d = 0; d < D; d++) buff[d] -= X[m * D + d];
00688                 DD[m] = .0;
00689                 for(int d = 0; d < D; d++) DD[m] += buff[d] * buff[d];
00690             }
00691             
00692             // Initialize some variables
00693             bool found = false;
00694             double beta = 1.0;
00695             double min_beta = -DBL_MAX;
00696             double max_beta =  DBL_MAX;
00697             double tol = 1e-5;
00698             
00699             // Iterate until we found a good perplexity
00700             int iter = 0; double sum_P;
00701             while(!found && iter < 200) {
00702                 
00703                 // Compute Gaussian kernel row
00704                 for(int m = 0; m < N; m++) cur_P[m] = exp(-beta * DD[m]);
00705                 cur_P[n] = DBL_MIN;
00706                 
00707                 // Compute entropy of current row
00708                 sum_P = DBL_MIN;
00709                 for(int m = 0; m < N; m++) sum_P += cur_P[m];
00710                 double H = 0.0;
00711                 for(int m = 0; m < N; m++) H += beta * (DD[m] * cur_P[m]);
00712                 H = (H / sum_P) + log(sum_P);
00713                 
00714                 // Evaluate whether the entropy is within the tolerance level
00715                 double Hdiff = H - log(perplexity);
00716                 if(Hdiff < tol && -Hdiff < tol) {
00717                     found = true;
00718                 }
00719                 else {
00720                     if(Hdiff > 0) {
00721                         min_beta = beta;
00722                         if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
00723                             beta *= 2.0;
00724                         else
00725                             beta = (beta + max_beta) / 2.0;
00726                     }
00727                     else {
00728                         max_beta = beta;
00729                         if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
00730                             beta /= 2.0;
00731                         else
00732                             beta = (beta + min_beta) / 2.0;
00733                     }
00734                 }
00735                 
00736                 // Update iteration counter
00737                 iter++;
00738             }
00739             
00740             // Row-normalize and threshold current row of P
00741             for(int m = 0; m < N; m++) cur_P[m] /= sum_P;
00742             for(int m = 0; m < N; m++) {
00743                 if(cur_P[m] > threshold / (double) N) {
00744                     col_P[count] = m;
00745                     val_P[count] = cur_P[m];
00746                     count++;
00747                 }
00748             }
00749             row_P[n + 1] = count;
00750         }
00751         
00752         // Clean up memory
00753         free(DD);    DD = NULL;
00754         free(buff);  buff = NULL;
00755         free(cur_P); cur_P = NULL;
00756     }
00757 
00758     void computeSquaredEuclideanDistance(double* X, int N, int D, double* DD)
00759     {
00760         double* dataSums = (double*) calloc(N, sizeof(double));
00761         if(dataSums == NULL) { printf("Memory allocation failed!\n"); exit(1); }
00762         for(int n = 0; n < N; n++) {
00763             for(int d = 0; d < D; d++) {
00764                 dataSums[n] += (X[n * D + d] * X[n * D + d]);
00765             }
00766         }
00767         for(int n = 0; n < N; n++) {
00768             for(int m = 0; m < N; m++) {
00769                 DD[n * N + m] = dataSums[n] + dataSums[m];
00770             }
00771         }
00772         Eigen::Map<Eigen::MatrixXd> DD_map(DD,N,N);
00773         Eigen::Map<Eigen::MatrixXd> X_map(X,D,N);
00774         DD_map.noalias() = -2.0*X_map.transpose()*X_map;
00775 
00776         //cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans, N, N, D, -2.0, X, D, X, D, 1.0, DD, N);
00777         free(dataSums); dataSums = NULL;
00778     }
00779 
00780 };
00781 
00782 }
00783 
00784 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines