From bd10ce4c56b4c6a44e2bd145c9200163a2b51f0d Mon Sep 17 00:00:00 2001 From: linqiaozhi Date: Thu, 7 Feb 2019 16:05:39 -0500 Subject: [PATCH 1/9] Supports variable degree of freedom for heavier tailed t-SNE Co-authored-by: dkobak --- fast_tsne.R | 5 +- fast_tsne.py | 3 +- src/nbodyfft.cpp | 4 +- src/nbodyfft.h | 4 +- src/tsne.cpp | 248 +++++++++++++++++++++++++++++++++++++++++------ src/tsne.h | 16 +-- 6 files changed, 240 insertions(+), 40 deletions(-) diff --git a/fast_tsne.R b/fast_tsne.R index 6c196a1..478b9f1 100644 --- a/fast_tsne.R +++ b/fast_tsne.R @@ -19,7 +19,7 @@ fftRtsne <- function(X, data_path=NULL, result_path=NULL, load_affinities=NULL, fast_tsne_path=NULL, nthreads=0, perplexity_list = NULL, - get_costs = FALSE, ... ) { + get_costs = FALSE, df = 1.0,... ) { if (is.null(fast_tsne_path)) { if(.Platform$OS.type == "unix") { @@ -51,6 +51,7 @@ fftRtsne <- function(X, if (!(max_iter>0)) { stop("Incorrect number of iterations.")} if (!is.wholenumber(stop_early_exag_iter) || stop_early_exag_iter<0) { stop("stop_early_exag_iter should be a positive integer")} if (!is.numeric(exaggeration_factor)) { stop("exaggeration_factor should be numeric")} + if (!is.numeric(df)) { stop("df should be numeric")} if (!is.wholenumber(dims) || dims<=0) { stop("Incorrect dimensionality.")} if (search_k == -1) { if (perplexity>0) { @@ -124,8 +125,10 @@ fftRtsne <- function(X, tX = c(t(X)) writeBin( tX, f) writeBin( as.integer(rand_seed), f,size=4) + writeBin(as.numeric(df), f, size=8) writeBin( as.integer(load_affinities), f,size=4) if (! is.null(initialization)){ writeBin( c(t(initialization)), f) } + print(df) close(f) flag= system2(command=fast_tsne_path, args=c(data_path, result_path, nthreads)); diff --git a/fast_tsne.py b/fast_tsne.py index 0a6e11b..e513df4 100644 --- a/fast_tsne.py +++ b/fast_tsne.py @@ -21,7 +21,7 @@ def fast_tsne(X, theta=.5, perplexity=30, map_dims=2, max_iter=1000, search_k=None, start_late_exag_iter=-1, late_exag_coeff=-1, nterms=3, intervals_per_integer=1, min_num_intervals=50, seed=-1, initialization=None, load_affinities=None, - perplexity_list=None): + perplexity_list=None, df=1): # X should be a numpy array of 64-bit doubles X = np.array(X).astype(float) @@ -90,6 +90,7 @@ def fast_tsne(X, theta=.5, perplexity=30, map_dims=2, max_iter=1000, f.write(struct.pack('=i', min_num_intervals)) f.write(X.tobytes()) f.write(struct.pack('=i', seed)) + f.write(struct.pack('=d', df)) f.write(struct.pack('=i', load_affinities)) if initialization is not None: diff --git a/src/nbodyfft.cpp b/src/nbodyfft.cpp index 4c1d9ee..7a7f8d4 100644 --- a/src/nbodyfft.cpp +++ b/src/nbodyfft.cpp @@ -8,7 +8,7 @@ void precompute_2d(double x_max, double x_min, double y_max, double y_min, int n_boxes, int n_interpolation_points, kernel_type_2d kernel, double *box_lower_bounds, double *box_upper_bounds, double *y_tilde_spacings, - double *y_tilde, double *x_tilde, complex *fft_kernel_tilde) { + double *y_tilde, double *x_tilde, complex *fft_kernel_tilde, double df ) { /* * Set up the boxes */ @@ -52,7 +52,7 @@ void precompute_2d(double x_max, double x_min, double y_max, double y_min, int n auto *kernel_tilde = new double[n_fft_coeffs * n_fft_coeffs](); for (int i = 0; i < n_interpolation_points_1d; i++) { for (int j = 0; j < n_interpolation_points_1d; j++) { - double tmp = kernel(y_tilde[0], x_tilde[0], y_tilde[i], x_tilde[j]); + double tmp = kernel(y_tilde[0], x_tilde[0], y_tilde[i], x_tilde[j],df ); kernel_tilde[(n_interpolation_points_1d + i) * n_fft_coeffs + (n_interpolation_points_1d + j)] = tmp; kernel_tilde[(n_interpolation_points_1d - i) * n_fft_coeffs + (n_interpolation_points_1d + j)] = tmp; kernel_tilde[(n_interpolation_points_1d + i) * n_fft_coeffs + (n_interpolation_points_1d - j)] = tmp; diff --git a/src/nbodyfft.h b/src/nbodyfft.h index 9127349..f46c4b3 100644 --- a/src/nbodyfft.h +++ b/src/nbodyfft.h @@ -12,11 +12,11 @@ using namespace std; typedef double (*kernel_type)(double, double); -typedef double (*kernel_type_2d)(double, double, double, double); +typedef double (*kernel_type_2d)(double, double, double, double, double); void precompute_2d(double x_max, double x_min, double y_max, double y_min, int n_boxes, int n_interpolation_points, kernel_type_2d kernel, double *box_lower_bounds, double *box_upper_bounds, double *y_tilde_spacings, - double *y_tilde, double *x_tilde, complex *fft_kernel_tilde); + double *y_tilde, double *x_tilde, complex *fft_kernel_tilde, double df); void n_body_fft_2d(int N, int n_terms, double *xs, double *ys, double *chargesQij, int n_boxes, int n_interpolation_points, double *box_lower_bounds, double *box_upper_bounds, diff --git a/src/tsne.cpp b/src/tsne.cpp index 1e5669a..374d076 100644 --- a/src/tsne.cpp +++ b/src/tsne.cpp @@ -71,10 +71,18 @@ double squared_cauchy(double x, double y) { } -double squared_cauchy_2d(double x1, double x2, double y1, double y2) { +double squared_cauchy_2d(double x1, double x2, double y1, double y2,double df) { return pow(1.0 + pow(x1 - y1, 2) + pow(x2 - y2, 2), -2); } +double general_kernel_2d(double x1, double x2, double y1, double y2, double df) { + return pow(1.0 + ((x1 - y1)*(x1-y1) + (x2 - y2)*(x2-y2))/df, -(df)); +} + +double squared_general_kernel_2d(double x1, double x2, double y1, double y2, double df) { + return pow(1.0 + ((x1 - y1)*(x1-y1) + (x2 - y2)*(x2-y2))/df, -(df+1.0)); +} + using namespace std; @@ -103,7 +111,7 @@ int TSNE::run(double *X, int N, int D, double *Y, int no_dims, double perplexity int nbody_algorithm, int knn_algo, double early_exag_coeff, double *costs, bool no_momentum_during_exag, int start_late_exag_iter, double late_exag_coeff, int n_trees, int search_k, int nterms, double intervals_per_integer, int min_num_intervals, unsigned int nthreads, - int load_affinities, int perplexity_list_length, double *perplexity_list) { + int load_affinities, int perplexity_list_length, double *perplexity_list, double df) { // Some logging messages if (N - 1 < 3 * perplexity) { @@ -420,7 +428,7 @@ int TSNE::run(double *X, int N, int D, double *Y, int no_dims, double perplexity if (exact) { // Compute the exact gradient using full P matrix - computeExactGradient(P, Y, N, no_dims, dY); + computeExactGradient(P, Y, N, no_dims, dY,df); } else { if (nbody_algorithm == 2) { // Use FFT accelerated interpolation based negative gradients @@ -428,8 +436,14 @@ int TSNE::run(double *X, int N, int D, double *Y, int no_dims, double perplexity computeFftGradientOneD(P, row_P, col_P, val_P, Y, N, no_dims, dY, nterms, intervals_per_integer, min_num_intervals, nthreads); } else { - computeFftGradient(P, row_P, col_P, val_P, Y, N, no_dims, dY, nterms, intervals_per_integer, - min_num_intervals, nthreads); + if (df ==1.0) { + computeFftGradient(P, row_P, col_P, val_P, Y, N, no_dims, dY, nterms, intervals_per_integer, + min_num_intervals, nthreads); + }else { + computeFftGradientVariableDf(P, row_P, col_P, val_P, Y, N, no_dims, dY, nterms, intervals_per_integer, + min_num_intervals, nthreads,df ); + + } } } else if (nbody_algorithm == 1) { // Otherwise, compute the negative gradient using the Barnes-Hut approximation @@ -441,7 +455,7 @@ int TSNE::run(double *X, int N, int D, double *Y, int no_dims, double perplexity computeGradient(P, row_P, col_P, val_P, Y, N, no_dims, dY, theta); computeFftGradient(P, row_P, col_P, val_P, Y, N, no_dims, dY, nterms, intervals_per_integer, min_num_intervals, nthreads); - computeExactGradientTest(Y, N, no_dims); + computeExactGradientTest(Y, N, no_dims,df); } // We can turn off momentum/gains until after the early exaggeration phase is completed @@ -487,10 +501,10 @@ int TSNE::run(double *X, int N, int D, double *Y, int no_dims, double perplexity START_TIME; double C = .0; if (exact) { - C = evaluateError(P, Y, N, no_dims); + C = evaluateError(P, Y, N, no_dims,df); }else{ if (nbody_algorithm == 2) { - C = evaluateErrorFft(row_P, col_P, val_P, Y, N, no_dims,nthreads); + C = evaluateErrorFft(row_P, col_P, val_P, Y, N, no_dims,nthreads,df); }else { C = evaluateError(row_P, col_P, val_P, Y, N, no_dims,theta, nthreads); } @@ -676,6 +690,173 @@ void TSNE::computeFftGradientOneD(double *P, unsigned int *inp_row_P, unsigned i delete[] neg_f; } +// Compute the gradient of the t-SNE cost function using the FFT interpolation +// based approximation, with variable degree of freedom df +void TSNE::computeFftGradientVariableDf(double *P, unsigned int *inp_row_P, unsigned int *inp_col_P, double *inp_val_P, double *Y, + int N, int D, double *dC, int n_interpolation_points, double intervals_per_integer, + int min_num_intervals, unsigned int nthreads, double df) { + + + // Zero out the gradient + for (int i = 0; i < N * D; i++) dC[i] = 0.0; + + // For convenience, split the x and y coordinate values + auto *xs = new double[N]; + auto *ys = new double[N]; + + double min_coord = INFINITY; + double max_coord = -INFINITY; + // Find the min/max values of the x and y coordinates + for (unsigned long i = 0; i < N; i++) { + xs[i] = Y[i * 2 + 0]; + ys[i] = Y[i * 2 + 1]; + if (xs[i] > max_coord) max_coord = xs[i]; + else if (xs[i] < min_coord) min_coord = xs[i]; + if (ys[i] > max_coord) max_coord = ys[i]; + else if (ys[i] < min_coord) min_coord = ys[i]; + } + // Compute the number of boxes in a single dimension and the total number of boxes in 2d + auto n_boxes_per_dim = static_cast(fmax(min_num_intervals, (max_coord - min_coord) / intervals_per_integer)); + + + //printf("min_coord: %lf, max_coord: %lf, n_boxes_per_dim: %d, (max_coord - min_coord) / intervals_per_integer) %d\n", min_coord, max_coord, n_boxes_per_dim, static_cast( (max_coord - min_coord) / intervals_per_integer)); + // FFTW works faster on numbers that can be written as 2^a 3^b 5^c 7^d + // 11^e 13^f, where e+f is either 0 or 1, and the other exponents are + // arbitrary + int allowed_n_boxes_per_dim[20] = {25,36, 50, 55, 60, 65, 70, 75, 80, 85, 90, 96, 100, 110, 120, 130, 140,150, 175, 200}; + if ( n_boxes_per_dim < allowed_n_boxes_per_dim[19] ) { + //Round up to nearest grid point + int chosen_i; + for (chosen_i =0; allowed_n_boxes_per_dim[chosen_i]< n_boxes_per_dim; chosen_i++); + n_boxes_per_dim = allowed_n_boxes_per_dim[chosen_i]; + } + + //printf(" n_boxes_per_dim: %d\n", n_boxes_per_dim ); + // The number of "charges" or s+2 sums i.e. number of kernel sums + int squared_n_terms = 3; + auto *SquaredChargesQij = new double[N * squared_n_terms]; + auto *SquaredPotentialsQij = new double[N * squared_n_terms](); + + // Prepare the terms that we'll use to compute the sum i.e. the repulsive forces + for (unsigned long j = 0; j < N; j++) { + SquaredChargesQij[j * squared_n_terms + 0] = xs[j]; + SquaredChargesQij[j * squared_n_terms + 1] = ys[j]; + SquaredChargesQij[j * squared_n_terms + 2] = 1; + } + + // Compute the number of boxes in a single dimension and the total number of boxes in 2d + int n_boxes = n_boxes_per_dim * n_boxes_per_dim; + + auto *box_lower_bounds = new double[2 * n_boxes]; + auto *box_upper_bounds = new double[2 * n_boxes]; + auto *y_tilde_spacings = new double[n_interpolation_points]; + int n_interpolation_points_1d = n_interpolation_points * n_boxes_per_dim; + auto *x_tilde = new double[n_interpolation_points_1d](); + auto *y_tilde = new double[n_interpolation_points_1d](); + auto *fft_kernel_tilde = new complex[2 * n_interpolation_points_1d * 2 * n_interpolation_points_1d]; + + INITIALIZE_TIME; + START_TIME; + precompute_2d(max_coord, min_coord, max_coord, min_coord, n_boxes_per_dim, n_interpolation_points, + &squared_general_kernel_2d, + box_lower_bounds, box_upper_bounds, y_tilde_spacings, x_tilde, y_tilde, fft_kernel_tilde, df); + n_body_fft_2d(N, squared_n_terms, xs, ys, SquaredChargesQij, n_boxes_per_dim, n_interpolation_points, box_lower_bounds, + box_upper_bounds, y_tilde_spacings, fft_kernel_tilde, SquaredPotentialsQij, nthreads); + + int not_squared_n_terms = 1; + auto *NotSquaredChargesQij = new double[N * not_squared_n_terms]; + auto *NotSquaredPotentialsQij = new double[N * not_squared_n_terms](); + + // Prepare the terms that we'll use to compute the sum i.e. the repulsive forces + for (unsigned long j = 0; j < N; j++) { + NotSquaredChargesQij[j * not_squared_n_terms + 0] = 1; + } + + precompute_2d(max_coord, min_coord, max_coord, min_coord, n_boxes_per_dim, n_interpolation_points, + &general_kernel_2d, + box_lower_bounds, box_upper_bounds, y_tilde_spacings, x_tilde, y_tilde, fft_kernel_tilde,df); + n_body_fft_2d(N, not_squared_n_terms, xs, ys, NotSquaredChargesQij, n_boxes_per_dim, n_interpolation_points, box_lower_bounds, + box_upper_bounds, y_tilde_spacings, fft_kernel_tilde, NotSquaredPotentialsQij, nthreads); + + + + + // Compute the normalization constant Z or sum of q_{ij}. + double sum_Q = 0; + for (unsigned long i = 0; i < N; i++) { + double h1 = NotSquaredPotentialsQij[i * not_squared_n_terms+ 0]; + sum_Q += h1; + } + sum_Q -= N; + + // Now, figure out the Gaussian component of the gradient. This corresponds to the "attraction" term of the + // gradient. It was calculated using a fast KNN approach, so here we just use the results that were passed to this + // function + unsigned int ind2 = 0; + double *pos_f = new double[N * 2]; + END_TIME("Total Interpolation"); + START_TIME; + // Loop over all edges in the graph + PARALLEL_FOR(nthreads, N, { + double dim1 = 0; + double dim2 = 0; + + for (unsigned int i = inp_row_P[loop_i]; i < inp_row_P[loop_i + 1]; i++) { + // Compute pairwise distance and Q-value + unsigned int ind3 = inp_col_P[i]; + double d_ij = (xs[loop_i] - xs[ind3]) * (xs[loop_i] - xs[ind3]) + (ys[loop_i] - ys[ind3]) * (ys[loop_i] - ys[ind3]); + double q_ij = 1 / (1 + d_ij/df); + + dim1 += inp_val_P[i] * q_ij * (xs[loop_i] - xs[ind3]); + dim2 += inp_val_P[i] * q_ij * (ys[loop_i] - ys[ind3]); + } + pos_f[loop_i * 2 + 0] = dim1; + pos_f[loop_i * 2 + 1] = dim2; + }); + + // Make the negative term, or F_rep in the equation 3 of the paper + END_TIME("Attractive Forces"); + + double *neg_f = new double[N * 2]; + for (unsigned int i = 0; i < N; i++) { + double h2 = SquaredPotentialsQij[i * squared_n_terms]; + double h3 = SquaredPotentialsQij[i * squared_n_terms + 1]; + double h4 = SquaredPotentialsQij[i * squared_n_terms + 2]; + neg_f[i * 2 + 0] = ( xs[i] *h4 - h2 ) / sum_Q; + neg_f[i * 2 + 1] = (ys[i] *h4 - h3 ) / sum_Q; + + dC[i * 2 + 0] = (pos_f[i * 2] - neg_f[i * 2]); + dC[i * 2 + 1] = (pos_f[i * 2 + 1] - neg_f[i * 2 + 1]); + + + } + + this->current_sum_Q = sum_Q; + +/* FILE *fp = nullptr; + char buffer[500]; + sprintf(buffer, "temp/fft_gradient%d.txt", itTest); + fp = fopen(buffer, "w"); // Open file for writing + for (int i = 0; i < N; i++) { + fprintf(fp, "%d,%.12e,%.12e\n", i, neg_f[i * 2] , neg_f[i * 2 + 1]); + } + fclose(fp);*/ + + delete[] pos_f; + delete[] neg_f; + delete[] SquaredPotentialsQij; + delete[] NotSquaredPotentialsQij; + delete[] SquaredChargesQij; + delete[] NotSquaredChargesQij; + delete[] xs; + delete[] ys; + delete[] box_lower_bounds; + delete[] box_upper_bounds; + delete[] y_tilde_spacings; + delete[] y_tilde; + delete[] x_tilde; + delete[] fft_kernel_tilde; +} // Compute the gradient of the t-SNE cost function using the FFT interpolation based approximation void TSNE::computeFftGradient(double *P, unsigned int *inp_row_P, unsigned int *inp_col_P, double *inp_val_P, double *Y, @@ -745,7 +926,7 @@ void TSNE::computeFftGradient(double *P, unsigned int *inp_row_P, unsigned int * START_TIME; precompute_2d(max_coord, min_coord, max_coord, min_coord, n_boxes_per_dim, n_interpolation_points, &squared_cauchy_2d, - box_lower_bounds, box_upper_bounds, y_tilde_spacings, x_tilde, y_tilde, fft_kernel_tilde); + box_lower_bounds, box_upper_bounds, y_tilde_spacings, x_tilde, y_tilde, fft_kernel_tilde,1.0); n_body_fft_2d(N, n_terms, xs, ys, chargesQij, n_boxes_per_dim, n_interpolation_points, box_lower_bounds, box_upper_bounds, y_tilde_spacings, fft_kernel_tilde, potentialsQij,nthreads); @@ -823,8 +1004,8 @@ void TSNE::computeFftGradient(double *P, unsigned int *inp_row_P, unsigned int * } -void TSNE::computeExactGradientTest(double *Y, int N, int D) { - // Compute the squared Euclidean distance matrix +void TSNE::computeExactGradientTest(double *Y, int N, int D, double df ) { + // Compute the squared Euclidean distance matrix double *DD = (double *) malloc(N * N * sizeof(double)); if (DD == NULL) { printf("Memory allocation failed!\n"); @@ -843,7 +1024,7 @@ void TSNE::computeExactGradientTest(double *Y, int N, int D) { for (int n = 0; n < N; n++) { for (int m = 0; m < N; m++) { if (n != m) { - Q[nN + m] = 1 / (1 + DD[nN + m]); + Q[nN + m] = 1.0 / pow(1.0 + DD[nN + m]/(double)df, (df)); sum_Q += Q[nN + m]; } } @@ -859,27 +1040,31 @@ void TSNE::computeExactGradientTest(double *Y, int N, int D) { for (int n = 0; n < N; n++) { double testQij = 0; double testPos = 0; - double testNeg = 0; + double testNeg1 = 0; + double testNeg2 = 0; double testdC = 0; int mD = 0; for (int m = 0; m < N; m++) { if (n != m) { - testNeg += Q[nN + m] * Q[nN + m] * (Y[nD + 0] - Y[mD + 0]) / sum_Q; + testNeg1 += pow(Q[nN + m],(df +1.0)/df) * (Y[nD + 0] - Y[mD + 0]) / sum_Q; + testNeg2 += pow(Q[nN + m],(df +1.0)/df) * (Y[nD + 1] - Y[mD + 1]) / sum_Q; } mD += D; } - fprintf(fp, "%d, %.12e\n", n, testNeg); + fprintf(fp, "%d, %.12e, %.12e\n", n, testNeg1,testNeg2); + nN += N; nD += D; } fclose(fp); free(DD); free(Q); + } // Compute the exact gradient of the t-SNE cost function -void TSNE::computeExactGradient(double *P, double *Y, int N, int D, double *dC) { +void TSNE::computeExactGradient(double *P, double *Y, int N, int D, double *dC, double df) { // Make sure the current gradient contains zeros for (int i = 0; i < N * D; i++) dC[i] = 0.0; @@ -897,7 +1082,7 @@ void TSNE::computeExactGradient(double *P, double *Y, int N, int D, double *dC) for (int n = 0; n < N; n++) { for (int m = 0; m < N; m++) { if (n != m) { - Q[nN + m] = 1 / (1 + DD[nN + m]); + Q[nN + m] = 1.0 / pow(1.0 + DD[nN + m]/(double)df, df); sum_Q += Q[nN + m]; } } @@ -911,7 +1096,7 @@ void TSNE::computeExactGradient(double *P, double *Y, int N, int D, double *dC) int mD = 0; for (int m = 0; m < N; m++) { if (n != m) { - double mult = (P[nN + m] - (Q[nN + m] / sum_Q)) * Q[nN + m]; + double mult = (P[nN + m] - (Q[nN + m] / sum_Q)) * pow(Q[nN + m], 1.0/df); for (int d = 0; d < D; d++) { dC[nD + d] += (Y[nD + d] - Y[mD + d]) * mult; } @@ -927,7 +1112,7 @@ void TSNE::computeExactGradient(double *P, double *Y, int N, int D, double *dC) // Evaluate t-SNE cost function (exactly) -double TSNE::evaluateError(double *P, double *Y, int N, int D) { +double TSNE::evaluateError(double *P, double *Y, int N, int D, double df) { // Compute the squared Euclidean distance matrix double *DD = (double *) malloc(N * N * sizeof(double)); double *Q = (double *) malloc(N * N * sizeof(double)); @@ -943,13 +1128,17 @@ double TSNE::evaluateError(double *P, double *Y, int N, int D) { for (int n = 0; n < N; n++) { for (int m = 0; m < N; m++) { if (n != m) { - Q[nN + m] = 1 / (1 + DD[nN + m]); + Q[nN + m] = 1.0 / pow(1.0 + DD[nN + m]/(double)df, df); sum_Q += Q[nN + m]; } else Q[nN + m] = DBL_MIN; } nN += N; } + //printf("sum_Q: %e", sum_Q); for (int i = 0; i < N * N; i++) Q[i] /= sum_Q; + // for (int i = 0; i < N; i++) printf("Q[%d]: %e\n", i, Q[i]); + +//printf("Q[N*N/2+1]: %e, Q[N*N-1]: %e\n", Q[N*N/2+1], Q[N*N/2+2]); // Sum t-SNE error double C = .0; @@ -964,7 +1153,7 @@ double TSNE::evaluateError(double *P, double *Y, int N, int D) { } // Evaluate t-SNE cost function (approximately) using FFT -double TSNE::evaluateErrorFft(unsigned int *row_P, unsigned int *col_P, double *val_P, double *Y, int N, int D,unsigned int nthreads) { +double TSNE::evaluateErrorFft(unsigned int *row_P, unsigned int *col_P, double *val_P, double *Y, int N, int D,unsigned int nthreads, double df) { // Get estimate of normalization term double sum_Q = this->current_sum_Q; @@ -981,7 +1170,7 @@ double TSNE::evaluateErrorFft(unsigned int *row_P, unsigned int *col_P, double * for (int d = 0; d < D; d++) buff[d] = Y[ind1 + d]; for (int d = 0; d < D; d++) buff[d] -= Y[ind2 + d]; for (int d = 0; d < D; d++) Q += buff[d] * buff[d]; - Q = (1.0 / (1.0 + Q)) / sum_Q; + Q = pow(1.0 / (1.0 + Q/df), df) / sum_Q; temp += val_P[i] * log((val_P[i] + FLT_MIN) / (Q + FLT_MIN)); } C += temp; @@ -1539,7 +1728,7 @@ bool TSNE::load_data(const char *data_path, double **data, double **Y, int *n, int *start_late_exag_iter, double *late_exag_coeff, int *nterms, double *intervals_per_integer, int *min_num_intervals, bool *skip_random_init, int *load_affinities, - int *perplexity_list_length, double **perplexity_list) { + int *perplexity_list_length, double **perplexity_list, double * df) { FILE *h; if((h = fopen(data_path, "r+b")) == NULL) { @@ -1598,6 +1787,9 @@ bool TSNE::load_data(const char *data_path, double **data, double **Y, int *n, if(!feof(h)) { result = fread(rand_seed, sizeof(int), 1, h); // random seed } + if(!feof(h)) { + result = fread(df, sizeof(double),1,h); // Number of degrees of freedom of the kernel + } if(!feof(h)) { result = fread(load_affinities, sizeof(int), 1, h); // to load or to save affinities } @@ -1616,7 +1808,6 @@ bool TSNE::load_data(const char *data_path, double **data, double **Y, int *n, } else{ *skip_random_init = false; } - fclose(h); printf("Read the following parameters:\n\t n %d by d %d dataset, theta %lf,\n" "\t perplexity %lf, no_dims %d, max_iter %d,\n" @@ -1626,14 +1817,14 @@ bool TSNE::load_data(const char *data_path, double **data, double **Y, int *n, "\t knn_algo %d, early_exag_coeff %lf,\n" "\t no_momentum_during_exag %d, n_trees %d, search_k %d,\n" "\t start_late_exag_iter %d, late_exag_coeff %lf\n" - "\t nterms %d, interval_per_integer %lf, min_num_intervals %d\n", + "\t nterms %d, interval_per_integer %lf, min_num_intervals %d, t-dist df %lf\n", *n, *d, *theta, *perplexity, *no_dims, *max_iter,*stop_lying_iter, *mom_switch_iter, *momentum, *final_momentum, *learning_rate, *K, *sigma, *nbody_algo, *knn_algo, *early_exag_coeff, *no_momentum_during_exag, *n_trees, *search_k, *start_late_exag_iter, *late_exag_coeff, - *nterms, *intervals_per_integer, *min_num_intervals); + *nterms, *intervals_per_integer, *min_num_intervals, *df); printf("Read the %i x %i data matrix successfully. X[0,0] = %lf\n", *n, *d, *data[0]); @@ -1695,6 +1886,7 @@ int main(int argc, char *argv[]) { double *perplexity_list; int perplexity_list_length; + double df; data_path = "data.dat"; result_path = "result.dat"; @@ -1724,7 +1916,7 @@ int main(int argc, char *argv[]) { &n_trees, &search_k, &start_late_exag_iter, &late_exag_coeff, &nterms, &intervals_per_integer, &min_num_intervals, &skip_random_init, &load_affinities, - &perplexity_list_length, &perplexity_list)) { + &perplexity_list_length, &perplexity_list, &df)) { bool no_momentum_during_exag_bool = true; if (no_momentum_during_exag == 0) no_momentum_during_exag_bool = false; @@ -1737,7 +1929,7 @@ int main(int argc, char *argv[]) { stop_lying_iter, mom_switch_iter, momentum, final_momentum, learning_rate, K, sigma, nbody_algo, knn_algo, early_exag_coeff, costs, no_momentum_during_exag_bool, start_late_exag_iter, late_exag_coeff, n_trees,search_k, nterms, intervals_per_integer, min_num_intervals, nthreads, load_affinities, - perplexity_list_length, perplexity_list); + perplexity_list_length, perplexity_list, df); if (error_code < 0) { exit(error_code); diff --git a/src/tsne.h b/src/tsne.h index a24f686..eb892a4 100644 --- a/src/tsne.h +++ b/src/tsne.h @@ -45,7 +45,7 @@ class TSNE { int nbody_algorithm, int knn_algo, double early_exag_coeff, double *costs, bool no_momentum_during_exag, int start_late_exag_iter, double late_exag_coeff, int n_trees, int search_k, int nterms, double intervals_per_integer, int min_num_intervals, unsigned int nthreads, int load_affinities, - int perplexity_list_length, double *perplexity_list); + int perplexity_list_length, double *perplexity_list, double df ); bool load_data(const char *data_path, double **data, double **Y, int *n, int *d, int *no_dims, double *theta, double *perplexity, int *rand_seed, int *max_iter, int *stop_lying_iter, @@ -53,7 +53,7 @@ class TSNE { int *nbody_algo, int *knn_algo, double* early_exag_coeff, int *no_momentum_during_exag, int *n_trees, int *search_k, int *start_late_exag_iter, double *late_exag_coeff, int *nterms, double *intervals_per_integer, int *min_num_intervals, bool *skip_random_init, int *load_affinities, - int *perplexity_list_length, double **perplexity_list); + int *perplexity_list_length, double **perplexity_list,double *df); void save_data(const char *result_path, double *data, double *costs, int n, int d, int max_iter); @@ -64,6 +64,10 @@ class TSNE { void computeGradient(double *P, unsigned int *inp_row_P, unsigned int *inp_col_P, double *inp_val_P, double *Y, int N, int D, double *dC, double theta); + void computeFftGradientVariableDf(double *P, unsigned int *inp_row_P, unsigned int *inp_col_P, double *inp_val_P, double *Y, + int N, int D, double *dC, int n_interpolation_points, double intervals_per_integer, + int min_num_intervals, unsigned int nthreads, double df); + void computeFftGradient(double *P, unsigned int *inp_row_P, unsigned int *inp_col_P, double *inp_val_P, double *Y, int N, int D, double *dC, int n_interpolation_points, double intervals_per_integer, int min_num_intervals, unsigned int nthreads); @@ -72,16 +76,16 @@ class TSNE { double *Y, int N, int D, double *dC, int n_interpolation_points, double intervals_per_integer, int min_num_intervals, unsigned int nthreads); - void computeExactGradient(double *P, double *Y, int N, int D, double *dC); + void computeExactGradient(double *P, double *Y, int N, int D, double *dC, double df); - void computeExactGradientTest(double *Y, int N, int D); + void computeExactGradientTest(double *Y, int N, int D, double df); - double evaluateError(double *P, double *Y, int N, int D); + double evaluateError(double *P, double *Y, int N, int D, double df); double evaluateError(unsigned int *row_P, unsigned int *col_P, double *val_P, double *Y, int N, int D, double theta, unsigned int nthreads); - double evaluateErrorFft(unsigned int *row_P, unsigned int *col_P, double *val_P, double *Y, int N, int D, unsigned int nthreads); + double evaluateErrorFft(unsigned int *row_P, unsigned int *col_P, double *val_P, double *Y, int N, int D, unsigned int nthreads, double df); void zeroMean(double *X, int N, int D); double distances2similarities(double *D, double *P, int N, int n, double perplexity, double sigma, bool ifSquared); From 01e181e787eae4c32b25a9dea531075227258333 Mon Sep 17 00:00:00 2001 From: linqiaozhi Date: Thu, 7 Feb 2019 19:18:15 -0500 Subject: [PATCH 2/9] Three commits below from kernel-independent2 -Fixed a mistake in the R wrapper support of returning the losses -Added return_loss support. Added nthreads support. -Removed duplicate pows Co-authored-by: dkobak --- fast_tsne.R | 1 + fast_tsne.py | 24 ++++++++++++++++++------ src/tsne.cpp | 16 ++++++++++++---- 3 files changed, 31 insertions(+), 10 deletions(-) diff --git a/fast_tsne.R b/fast_tsne.R index 478b9f1..01c4072 100644 --- a/fast_tsne.R +++ b/fast_tsne.R @@ -141,6 +141,7 @@ fftRtsne <- function(X, Y <- readBin(f, numeric(), n=n*d); Y <- t(matrix(Y, nrow=d)); if (get_costs ) { + tmp <- readBin(f, integer(), n=1, size=4); costs <- readBin(f, numeric(), n=max_iter,size=8); Yout <- list( Y=Y, costs=costs); }else { diff --git a/fast_tsne.py b/fast_tsne.py index e513df4..fbc4c89 100644 --- a/fast_tsne.py +++ b/fast_tsne.py @@ -2,10 +2,11 @@ # # Usage example: # import sys; sys.path.append('../') -# from fast_tsne import fast_tsne -# import numpy as np +# from fast_tsne import fast_tsne +# import numpy as np # X = np.random.randn(1000, 50) # Z = fast_tsne(X, perplexity = 30) +# # Written by Dmitry Kobak @@ -21,7 +22,7 @@ def fast_tsne(X, theta=.5, perplexity=30, map_dims=2, max_iter=1000, search_k=None, start_late_exag_iter=-1, late_exag_coeff=-1, nterms=3, intervals_per_integer=1, min_num_intervals=50, seed=-1, initialization=None, load_affinities=None, - perplexity_list=None, df=1): + perplexity_list=None, df=1, return_loss=False, nthreads=0): # X should be a numpy array of 64-bit doubles X = np.array(X).astype(float) @@ -98,15 +99,26 @@ def fast_tsne(X, theta=.5, perplexity=30, map_dims=2, max_iter=1000, f.write(initialization.tobytes()) # run t-sne - subprocess.call(os.path.dirname(os.path.realpath(__file__)) + '/bin/fast_tsne') + subprocess.call([os.path.dirname(os.path.realpath(__file__)) + + '/bin/fast_tsne', 'data.dat', 'result.dat', '{}'.format(nthreads)]) # read data file with open(os.getcwd()+'/result.dat', 'rb') as f: n, = struct.unpack('=i', f.read(4)) md, = struct.unpack('=i', f.read(4)) sz = struct.calcsize('=d') - buf = f.read() + buf = f.read(sz*n*md) x_tsne = [struct.unpack_from('=d', buf, sz*offset) for offset in range(n*md)] x_tsne = np.array(x_tsne).reshape((n,md)) + _, = struct.unpack('=i', f.read(4)) + buf = f.read(sz*max_iter) + loss = [struct.unpack_from('=d', buf, sz*offset) for offset in range(max_iter)] + loss = np.array(loss).squeeze() + loss[np.arange(1,max_iter+1)%50>0] = np.nan + + if return_loss: + return (x_tsne, loss) + else: + return x_tsne + - return x_tsne diff --git a/src/tsne.cpp b/src/tsne.cpp index 374d076..72825d9 100644 --- a/src/tsne.cpp +++ b/src/tsne.cpp @@ -1077,13 +1077,18 @@ void TSNE::computeExactGradient(double *P, double *Y, int N, int D, double *dC, auto *Q = (double *) malloc(N * N * sizeof(double)); if (Q == nullptr) throw std::bad_alloc(); + auto *Qpow = (double *) malloc(N * N * sizeof(double)); + if (Qpow == nullptr) throw std::bad_alloc(); + double sum_Q = .0; int nN = 0; for (int n = 0; n < N; n++) { for (int m = 0; m < N; m++) { if (n != m) { - Q[nN + m] = 1.0 / pow(1.0 + DD[nN + m]/(double)df, df); - sum_Q += Q[nN + m]; + //Q[nN + m] = 1.0 / pow(1.0 + DD[nN + m]/(double)df, df); + Q[nN + m] = 1.0 / (1.0 + DD[nN + m]/(double)df); + Qpow[nN + m] = pow(Q[nN + m], df); + sum_Q += Qpow[nN + m]; } } nN += N; @@ -1096,7 +1101,7 @@ void TSNE::computeExactGradient(double *P, double *Y, int N, int D, double *dC, int mD = 0; for (int m = 0; m < N; m++) { if (n != m) { - double mult = (P[nN + m] - (Q[nN + m] / sum_Q)) * pow(Q[nN + m], 1.0/df); + double mult = (P[nN + m] - (Qpow[nN + m] / sum_Q)) * (Q[nN + m]); for (int d = 0; d < D; d++) { dC[nD + d] += (Y[nD + d] - Y[mD + d]) * mult; } @@ -1107,6 +1112,7 @@ void TSNE::computeExactGradient(double *P, double *Y, int N, int D, double *dC, nD += D; } free(Q); + free(Qpow); free(DD); } @@ -1128,7 +1134,9 @@ double TSNE::evaluateError(double *P, double *Y, int N, int D, double df) { for (int n = 0; n < N; n++) { for (int m = 0; m < N; m++) { if (n != m) { - Q[nN + m] = 1.0 / pow(1.0 + DD[nN + m]/(double)df, df); + //Q[nN + m] = 1.0 / pow(1.0 + DD[nN + m]/(double)df, df); + Q[nN + m] = 1.0 / (1.0 + DD[nN + m]/(double)df); + Q[nN +m ] = pow(Q[nN +m ], df); sum_Q += Q[nN + m]; } else Q[nN + m] = DBL_MIN; } From 0f7e5c5f9a2c4e7cb5238002004d5f4363dcea81 Mon Sep 17 00:00:00 2001 From: linqiaozhi Date: Thu, 7 Feb 2019 22:48:41 -0500 Subject: [PATCH 3/9] Implemented version check in the MATLAB and R wrappers; fixed bug in MATLAB wrapper that was preventing it from running, introduced in #63 --- fast_tsne.R | 3 ++- fast_tsne.m | 17 +++++++++++++---- src/tsne.cpp | 25 ++++++++++++++++++------- 3 files changed, 33 insertions(+), 12 deletions(-) diff --git a/fast_tsne.R b/fast_tsne.R index 01c4072..58e46f8 100644 --- a/fast_tsne.R +++ b/fast_tsne.R @@ -20,6 +20,7 @@ fftRtsne <- function(X, load_affinities=NULL, fast_tsne_path=NULL, nthreads=0, perplexity_list = NULL, get_costs = FALSE, df = 1.0,... ) { + version_number = '1.1.0' if (is.null(fast_tsne_path)) { if(.Platform$OS.type == "unix") { @@ -131,7 +132,7 @@ fftRtsne <- function(X, print(df) close(f) - flag= system2(command=fast_tsne_path, args=c(data_path, result_path, nthreads)); + flag= system2(command=fast_tsne_path, args=c(version_number,data_path, result_path, nthreads)); if (flag != 0) { stop('tsne call failed'); } diff --git a/fast_tsne.m b/fast_tsne.m index b6090a5..0989f70 100644 --- a/fast_tsne.m +++ b/fast_tsne.m @@ -96,6 +96,8 @@ % CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING % IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY % OF SUCH DAMAGE. + + version_number = '1.1.0'; if (nargin == 1) opts.perplexity = 30; end @@ -275,6 +277,12 @@ else nthreads = opts.nthreads; end + + if (~isfield(opts, 'df')) + df = 1; + else + df = opts.df; + end X = double(X); @@ -283,7 +291,7 @@ % Compile t-SNE C code if(~exist(fullfile(tsne_path,'./fast_tsne'),'file') && isunix) - system(sprintf('g++ -std=c++11 -O3 src/sptree.cpp src/tsne.cpp src/nbodyfft.cpp -o bin/fast_tsne -pthread -lfftw3 -lm')); end + system(sprintf('g++ -std=c++11 -O3 src/sptree.cpp src/tsne.cpp src/nbodyfft.cpp -o bin/fast_tsne -pthread -lfftw3 -lm')); end % Compile t-SNE C code on Windows @@ -296,12 +304,12 @@ stop_lying_iter, K, sigma, nbody_algo, no_momentum_during_exag, knn_algo,... early_exag_coeff, n_trees, search_k, start_late_exag_iter, late_exag_coeff, rand_seed,... nterms, intervals_per_integer, min_num_intervals, initialization, load_affinities, ... - perplexity_list, mom_switch_iter, momentum, final_momentum, learning_rate); + perplexity_list, mom_switch_iter, momentum, final_momentum, learning_rate,df); disp('Data written'); tic %[flag, cmdout] = system(fullfile(tsne_path,'/fast_tsne'), '-echo'); - cmd = sprintf('%s data.dat result.dat %d',fullfile(tsne_path,'/fast_tsne'), nthreads); + cmd = sprintf('%s %s data.dat result.dat %d',fullfile(tsne_path,'/fast_tsne'), version_number, nthreads); [flag, cmdout] = system(cmd, '-echo'); if(flag~=0) error(cmdout); @@ -318,7 +326,7 @@ function write_data(filename, X, no_dims, theta, perplexity, max_iter,... stop_lying_iter, K, sigma, nbody_algo, no_momentum_during_exag, knn_algo,... early_exag_coeff, n_trees, search_k, start_late_exag_iter, late_exag_coeff, rand_seed,... nterms, intervals_per_integer, min_num_intervals, initialization, load_affinities, ... - perplexity_list, mom_switch_iter, momentum, final_momentum, learning_rate) + perplexity_list, mom_switch_iter, momentum, final_momentum, learning_rate,df) [n, d] = size(X); @@ -353,6 +361,7 @@ function write_data(filename, X, no_dims, theta, perplexity, max_iter,... fwrite(h, min_num_intervals, 'int'); fwrite(h, X', 'double'); fwrite(h, rand_seed, 'integer*4'); + fwrite(h, df, 'double'); fwrite(h, load_affinities, 'integer*4'); if ~isnan(initialization) fwrite(h, initialization', 'double'); diff --git a/src/tsne.cpp b/src/tsne.cpp index 72825d9..fd01ada 100644 --- a/src/tsne.cpp +++ b/src/tsne.cpp @@ -1871,13 +1871,14 @@ void TSNE::save_data(const char *result_path, double* data, double* costs, int n int main(int argc, char *argv[]) { - printf("=============== t-SNE v1.0.1 ===============\n"); + const char version_number[] = "1.1.0"; + printf("=============== t-SNE v%s ===============\n", version_number); // Define some variables int N, D, no_dims, max_iter, stop_lying_iter; int K, nbody_algo, knn_algo, no_momentum_during_exag; - int mom_switch_iter; - double momentum, final_momentum, learning_rate; + int mom_switch_iter; + double momentum, final_momentum, learning_rate; int n_trees, search_k, start_late_exag_iter; double sigma, early_exag_coeff, late_exag_coeff; double perplexity, theta, *data, *initial_data; @@ -1899,13 +1900,23 @@ int main(int argc, char *argv[]) { data_path = "data.dat"; result_path = "result.dat"; nthreads = 0; - if(argc >= 2) { - data_path = argv[1]; - } + if (argc >=2 ) { + if ( strcmp(argv[1],version_number)) { + std::cout<<"Wrapper passed wrong version number: "<< argv[1] <= 3) { - result_path = argv[2]; + data_path = argv[2]; } if(argc >= 4) { + result_path = argv[3]; + } + if(argc >= 5) { nthreads = (unsigned int)strtoul(argv[3], (char **)NULL, 10); } if (nthreads == 0) { From 48b15f71fcee79298d88f05a58d1a9499fa582b1 Mon Sep 17 00:00:00 2001 From: linqiaozhi Date: Thu, 7 Feb 2019 23:40:58 -0500 Subject: [PATCH 4/9] Updated readme to better describe features --- README.md | 30 +++++++++++++++++++----------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/README.md b/README.md index 42e808f..f1aa370 100644 --- a/README.md +++ b/README.md @@ -1,17 +1,25 @@ # FFT-accelerated Interpolation-based t-SNE (FIt-SNE) ## Introduction -t-Stochastic Neighborhood Embedding ([t-SNE](https://lvdmaaten.github.io/tsne/)) is a highly successful method for dimensionality reduction and visualization of high dimensional datasets. A popular [implementation](https://github.com/lvdmaaten/bhtsne) of t-SNE uses the Barnes-Hut algorithm to approximate the gradient at each iteration of gradient descent. We modified this implementation as follows: +t-Stochastic Neighborhood Embedding ([t-SNE](https://lvdmaaten.github.io/tsne/)) is a highly successful method for dimensionality reduction and visualization of high dimensional datasets. A popular [implementation](https://github.com/lvdmaaten/bhtsne) of t-SNE uses the Barnes-Hut algorithm to approximate the gradient at each iteration of gradient descent. We accelerated this implementation as follows: * Computation of the N-body Simulation: Instead of approximating the N-body simulation using Barnes-Hut, we interpolate onto an equispaced grid and use FFT to perform the convolution, dramatically reducing the time to compute the gradient at each iteration of gradient descent. See the [this](http://gauss.math.yale.edu/~gcl22/blog/numerics/low-rank/t-sne/2018/01/11/low-rank-kernels.html) post for some intuition on how it works. -* Computation of Input Similiarities: Instead of computing nearest neighbors using vantage-point trees, we approximate nearest neighbors using the [Annoy](https://github.com/spotify/annoy) library. The neighbor lookups are multithreaded to take advantage of machines with multiple cores. Using "near" neighbors as opposed to strictly "nearest" neighbors is faster, but also has a smoothing effect, which can be useful for embedding some datasets (see [Linderman et al. (2017)](https://arxiv.org/abs/1711.04712)). If subtle detail is required (e.g. in identifying small clusters), then use vantage-point trees (which is also multithreaded in this implementation). -* Early exaggeration: In [Linderman and Steinerberger (2017)](https://arxiv.org/abs/1706.02582), we showed that appropriately choosing the early exaggeration coefficient can lead to improved embedding of swissrolls and other synthetic datasets. -* Late exaggeration: Increasing the exaggeration coefficient late in the optimization process (e.g. after 800 of 1000 iterations) can improve separation of the clusters. +* Computation of Input Similarities: Instead of computing nearest neighbors using vantage-point trees, we approximate nearest neighbors using the [Annoy](https://github.com/spotify/annoy) library. The neighbor lookups are multithreaded to take advantage of machines with multiple cores. Using "near" neighbors as opposed to strictly "nearest" neighbors is faster, but also has a smoothing effect, which can be useful for embedding some datasets (see [Linderman et al. (2017)](https://arxiv.org/abs/1711.04712)). If subtle detail is required (e.g. in identifying small clusters), then use vantage-point trees (which is also multithreaded in this implementation). + Check out our [preprint](https://arxiv.org/abs/1712.09005) for more details and some benchmarks. -R, Matlab, and Python wrappers are `fast_tsne.R`, `fast_tsne.m`, and `fast_tsne.py` respectively. [Gioele La Manno](https://twitter.com/GioeleLaManno) implemented a Python (Cython) wrapper, which is available on PyPI [here](https://pypi.python.org/pypi/fitsne). +## Features +Additionally, this implementation includes the following features: +* Early exaggeration: In [Linderman and Steinerberger (2017)](https://arxiv.org/abs/1706.02582), we showed that appropriately choosing the early exaggeration coefficient can lead to improved embedding of swissrolls and other synthetic datasets. Early exaggeration is built into all t-SNE implementations; here we highlight its importance as a parameter. +* Late exaggeration: Increasing the exaggeration coefficient late in the optimization process can improve separation of the clusters. For example, [Kobak and Berens (2018)](https://www.biorxiv.org/content/10.1101/453449v1) suggest starting late exaggeration immediately following early exaggeration. +* Perplexity annealing: The perplexity parameter determines the width of the Gaussian kernel, such that small perplexity values uncover local structure while larger values reveal global structure. [Kobak and Berens (2018)](https://www.biorxiv.org/content/10.1101/453449v1) suggest using a range of perplexity values, where the large perplexity is used as initialization for t-SNE with small perplexity, resulting in a multi-scale embedding. +* PCA initialization: As suggested by [Kobak and Berens (2018)](https://www.biorxiv.org/content/10.1101/453449v1), initializing t-SNE with the first two principal components (scaled to have standard deviation 0.0001) results in an embedding which preserves the global structure more effectively than the default random normalization. +* Downsampling-based initialization: +* Variable degrees of freedom: [Kobak et al. (2019)]() show that decreasing the degree of freedom (df) of the t-distribution (resulting in heavier tails) reveals fine structure that is not visible in standard t-SNE. + ## Installation +R, Matlab, and Python wrappers are `fast_tsne.R`, `fast_tsne.m`, and `fast_tsne.py` respectively. Each of these wrappers can be used after installing FFTW and compiling the C++ code, as below. [Gioele La Manno](https://twitter.com/GioeleLaManno) implemented a Python (Cython) wrapper, which is available on PyPI [here](https://pypi.python.org/pypi/fitsne). **Note:** If you update to a new version of FIt-SNE using `git pull`, be sure to recompile. @@ -34,18 +42,18 @@ If you would like to compile it yourself see below. The code has been currently 2. Copy the binary file ( e.g. `x64/Debug/FItSNE.exe`) generated by the build process to the `bin/` folder 3. For Windows, we have added all dependencies, including the [FFTW library](http://www.fftw.org/), which is distributed under the GNU General Public License. For the binary to find the FFTW DLLs, you need to either add `src/winlibs/fftw/` to your PATH, or to copy the DLLs into the `bin/` directory. -As of this commit, only the R wrapper properly calls the Windows executable. -The Python and Matlab wrappers can be trivially changed to call it (just -changing `bin/fast_tsne` to `bin/FItSNE.exe` in the code), and will be changed -in future commits. +As of this commit, only the R wrapper properly calls the Windows executable. The Python and Matlab wrappers can be trivially changed to call it (just changing `bin/fast_tsne` to `bin/FItSNE.exe` in the code), and will be changed in future commits. Many thanks to [Josef Spidlen](https://github.com/jspidlen) for this Windows implementation! -## References -If you use our software, please cite: +## Acknowledgements and References +We are grateful for members of the community who have [contributed](https://github.com/KlugerLab/FIt-SNE/graphs/contributors) to improving FIt-SNE, especially [Dmitry Kobak](https://github.com/dkobak), [Pavlin Poličar](https://github.com/pavlin-policar), and [Josef Spidlen](https://github.com/jspidlen). + +If you use FIt-SNE, please cite: George C. Linderman, Manas Rachh, Jeremy G. Hoskins, Stefan Steinerberger, Yuval Kluger. (2017). Efficient Algorithms for t-distributed Stochastic Neighborhood Embedding. (2017) *arXiv:1712.09005* ([link](https://arxiv.org/abs/1712.09005)) Our implementation is derived from the Barnes-Hut implementation: Laurens van der Maaten (2014). Accelerating t-SNE using tree-based algorithms. Journal of Machine Learning Research, 15(1):3221–3245. ([link](https://dl.acm.org/citation.cfm?id=2627435.2697068)) + From cdf92d226d66e6315fa83ea7316f495cfeba8ede Mon Sep 17 00:00:00 2001 From: linqiaozhi Date: Thu, 7 Feb 2019 23:52:05 -0500 Subject: [PATCH 5/9] Documentation in R and MATLAB wrappers --- fast_tsne.R | 64 +++++++++++++++++++++++++++++++++++++++++++++++++++++ fast_tsne.m | 5 +++++ 2 files changed, 69 insertions(+) diff --git a/fast_tsne.R b/fast_tsne.R index 58e46f8..1bbbeba 100644 --- a/fast_tsne.R +++ b/fast_tsne.R @@ -3,6 +3,70 @@ FAST_TSNE_SCRIPT_DIR <<-getwd() cat(sprintf("FIt-SNE R wrapper loading.\nFIt-SNE root directory was set to %s\n", FAST_TSNE_SCRIPT_DIR)) +# Compute FIt-SNE of a dataset +# dims - dimensionality of the embedding. Default 2. +# perplexity - perplexity is used to determine the +# bandwidth of the Gaussian kernel in the input +# space. Default 30. +# theta - Set to 0 for exact. If non-zero, then will use either +# Barnes Hut or FIt-SNE based on nbody_algo. If Barnes Hut, then +# this determins the accuracy of BH approximation. +# Default 0.5. +# max_iter - Number of iterations of t-SNE to run. +# Default 1000. +# fft_not_bh - if theta is nonzero, this determins whether to +# use FIt-SNE or Barnes Hut approximation. Default is FIt-SNE. +# set to be True for FIt-SNE +# ann_not_vptree - use vp-trees (as in bhtsne) or approximate nearest neighbors (default). +# set to be True for approximate nearest neighbors +# exaggeration_factor - coefficient for early exaggeration +# (>1). Default 12. +# no_momentum_during_exag - Set to 0 to use momentum +# and other optimization tricks. 1 to do plain,vanilla +# gradient descent (useful for testing large exaggeration +# coefficients) +# stop_early_exag_iter - When to switch off early exaggeration. +# Default 250. +# start_late_exag_iter - When to start late +# exaggeration. set to -1 to not use late exaggeration +# Default -1. +# late_exag_coeff - Late exaggeration coefficient. +# Set to -1 to not use late exaggeration. +# Default -1 +# nterms - If using FIt-SNE, this is the number of +# interpolation points per sub-interval +# intervals_per_integer - See min_num_intervals +# min_num_intervals - Let maxloc = ceil(max(max(X))) +# and minloc = floor(min(min(X))). i.e. the points are in +# a [minloc]^no_dims by [maxloc]^no_dims interval/square. +# The number of intervals in each dimension is either +# min_num_intervals or ceil((maxloc - +# minloc)/intervals_per_integer), whichever is +# larger. min_num_intervals must be an integer >0, +# and intervals_per_integer must be >0. Default: +# min_num_intervals=50, intervals_per_integer = +# 1 +# +# sigma - Fixed sigma value to use when perplexity==-1 +# Default -1 (None) +# K - Number of nearest neighbours to get when using fixed sigma +# Default -30 (None) +# +# initialization - N x no_dims array to intialize the solution +# Default: None +# +# load_affinities - +# If 1, input similarities are loaded from a file and not computed +# If 2, input similarities are saved into a file. +# If 0, affinities are neither saved nor loaded +# +# perplexity_list - if perplexity==0 then perplexity combination will +# be used with values taken from perplexity_list. Default: NULL +# df - Degree of freedom of t-distribution, must be greater than 0. +# Values smaller than 1 correspond to heavier tails, which can often +# resolve substructure in the embedding. See Kobak et al. (2019) for +# details. Default is 1.0 +# fftRtsne <- function(X, dims=2, perplexity=30, theta=0.5, check_duplicates=TRUE, diff --git a/fast_tsne.m b/fast_tsne.m index 0989f70..01ab96d 100644 --- a/fast_tsne.m +++ b/fast_tsne.m @@ -63,6 +63,11 @@ % % opts.perplexity_list - if perplexity==0 then perplexity combination will % be used with values taken from perplexity_list. Default: [] +% opts.df - Degree of freedom of t-distribution, must be greater than 0. +% Values smaller than 1 correspond to heavier tails, which can often +% resolve substructure in the embedding. See Kobak et al. (2019) for +% details. Default is 1.0 + % Runs the C++ implementation of fast t-SNE using either the IFt-SNE From 1cc2fb4a4bf893e28764d5bb7c5c2e686c14283e Mon Sep 17 00:00:00 2001 From: Dmitry Kobak Date: Fri, 8 Feb 2019 09:38:35 +0100 Subject: [PATCH 6/9] Added documentation to the Python code. Added version control to the Python code. Added df!=1 to the Python example notebook --- examples/test.ipynb | 122 +++++++++++++++++++++++++++++++++----------- fast_tsne.R | 2 +- fast_tsne.py | 116 +++++++++++++++++++++++++++++++++++++++-- 3 files changed, 206 insertions(+), 34 deletions(-) diff --git a/examples/test.ipynb b/examples/test.ipynb index 6cc61ec..065f2fb 100644 --- a/examples/test.ipynb +++ b/examples/test.ipynb @@ -24,19 +24,23 @@ "import seaborn as sns; sns.set()\n", "\n", "# change the path!\n", - "import sys; sys.path.append('/home/dmitry/github/FIt-SNE')\n", + "import sys; sys.path.append('/home/dmitry/github/kerind/FIt-SNE')\n", "from fast_tsne import fast_tsne" ] }, { "cell_type": "code", "execution_count": 2, - "metadata": {}, + "metadata": { + "collapsed": false + }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ + "/home/dmitry/anaconda3/lib/python3.6/site-packages/h5py/__init__.py:34: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n", + " from ._conv import register_converters as _register_converters\n", "Using TensorFlow backend.\n" ] }, @@ -73,8 +77,10 @@ }, { "cell_type": "code", - "execution_count": 4, - "metadata": {}, + "execution_count": 3, + "metadata": { + "collapsed": false + }, "outputs": [ { "data": { @@ -158,7 +164,7 @@ " };\n", "\n", " this.imageObj.onunload = function() {\n", - " this.ws.close();\n", + " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", @@ -653,7 +659,7 @@ " // Register the callback with on_msg.\n", " comm.on_msg(function(msg) {\n", " //console.log('receiving', msg['content']['data'], msg);\n", - " // Pass the mpl event to the overriden (by mpl) onmessage function.\n", + " // Pass the mpl event to the overridden (by mpl) onmessage function.\n", " ws.onmessage(msg['content']['data'])\n", " });\n", " return ws;\n", @@ -805,9 +811,12 @@ " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", - " // select the cell after this one\n", - " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", - " IPython.notebook.select(index + 1);\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", " }\n", "}\n", "\n", @@ -856,7 +865,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -879,7 +888,9 @@ { "cell_type": "code", "execution_count": 4, - "metadata": {}, + "metadata": { + "collapsed": false + }, "outputs": [ { "data": { @@ -1673,14 +1684,22 @@ ], "source": [ "# Two classes are separated into two parts on the above figure. \n", - "# This can be fixed if one uses stronger/longer early exaggeration.\n", + "# This can be fixed if one uses stronger/longer early exaggeration (e.g. stop_early_exag_iter=500)\n", + "# or higher learning rate (e.g. learning rate=1000)\n", "\n", "# PCA initialization is not important for MNIST, but can be handy for other datasets\n", "\n", - "Z = fast_tsne(X50, perplexity=50, stop_early_exag_iter=500, initialization=PCAinit)\n", + "Z1 = fast_tsne(X50, perplexity=50, stop_early_exag_iter=500, initialization=PCAinit)\n", "\n", - "plt.figure(figsize=(5,5))\n", - "plt.scatter(Z[:,0], Z[:,1], c=col[y], s=1)\n", + "Z2 = fast_tsne(X50, perplexity=50, learning_rate=1000, initialization=PCAinit)\n", + "\n", + "plt.figure(figsize=(8,4))\n", + "plt.subplot(121)\n", + "plt.scatter(Z1[:,0], Z1[:,1], c=col[y], s=1)\n", + "plt.title('Early exaggeration for 500 iterations')\n", + "plt.subplot(122)\n", + "plt.scatter(Z2[:,0], Z2[:,1], c=col[y], s=1)\n", + "plt.title('Learning rate 1000')\n", "plt.tight_layout()" ] }, @@ -1716,7 +1735,9 @@ { "cell_type": "code", "execution_count": 6, - "metadata": {}, + "metadata": { + "collapsed": false + }, "outputs": [ { "data": { @@ -2536,7 +2557,9 @@ { "cell_type": "code", "execution_count": 7, - "metadata": {}, + "metadata": { + "collapsed": false + }, "outputs": [ { "data": { @@ -3356,7 +3379,9 @@ { "cell_type": "code", "execution_count": 6, - "metadata": {}, + "metadata": { + "collapsed": false + }, "outputs": [ { "data": { @@ -4176,7 +4201,9 @@ { "cell_type": "code", "execution_count": 11, - "metadata": {}, + "metadata": { + "collapsed": false + }, "outputs": [ { "data": { @@ -4995,7 +5022,9 @@ { "cell_type": "code", "execution_count": 9, - "metadata": {}, + "metadata": { + "collapsed": false + }, "outputs": [ { "data": { @@ -5788,18 +5817,18 @@ } ], "source": [ - "Z1 = fast_tsne(X50, perplexity=50, stop_early_exag_iter=500, initialization=PCAinit)\n", + "Z1 = fast_tsne(X50, perplexity=50, learning_rate=1000, initialization=PCAinit)\n", "\n", - "Z2 = fast_tsne(X50, perplexity=50, stop_early_exag_iter=500, initialization=PCAinit, \n", - " late_exag_coeff=4, start_late_exag_iter=500)\n", + "Z2 = fast_tsne(X50, perplexity=50, learning_rate=1000, initialization=PCAinit, \n", + " late_exag_coeff=4, start_late_exag_iter=250)\n", "\n", "plt.figure(figsize=(8,4))\n", "plt.subplot(121)\n", "plt.scatter(Z1[:,0], Z1[:,1], c=col[y], s=1)\n", - "plt.title('Early exagg=12 until iter=500')\n", + "plt.title('Early exagg=12 until iter=250')\n", "plt.subplot(122)\n", "plt.scatter(Z2[:,0], Z2[:,1], c=col[y], s=1)\n", - "plt.title('Early exagg=12 until iter=500\\nLate exagg=4 from iter=500')\n", + "plt.title('Early exagg=12 until iter=250\\nLate exagg=4 from iter=250')\n", "plt.tight_layout()" ] }, @@ -5813,7 +5842,9 @@ { "cell_type": "code", "execution_count": 12, - "metadata": {}, + "metadata": { + "collapsed": false + }, "outputs": [ { "data": { @@ -7442,7 +7473,9 @@ { "cell_type": "code", "execution_count": 21, - "metadata": {}, + "metadata": { + "collapsed": false + }, "outputs": [ { "data": { @@ -8261,7 +8294,9 @@ { "cell_type": "code", "execution_count": 10, - "metadata": {}, + "metadata": { + "collapsed": false + }, "outputs": [ { "data": { @@ -9069,6 +9104,35 @@ "plt.title('Perplexity = [10, 50, 100]')\n", "plt.tight_layout()" ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Making the kernel more heavy-tailed" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "Z1 = fast_tsne(X50, perplexity=50, learning_rate=1000, initialization=PCAinit)\n", + "\n", + "Z2 = fast_tsne(X50, perplexity=50, learning_rate=1000, initialization=PCAinit, df=.5)\n", + "\n", + "plt.figure(figsize=(8,4))\n", + "plt.subplot(121)\n", + "plt.scatter(Z1[:,0], Z1[:,1], c=col[y], s=1)\n", + "plt.title('1/(1+x^2) kernel')\n", + "plt.subplot(122)\n", + "plt.scatter(Z2[:,0], Z2[:,1], c=col[y], s=1)\n", + "plt.title('1/(1+x^2/0.5)^0.5 kernel')\n", + "plt.tight_layout()" + ] } ], "metadata": { @@ -9087,7 +9151,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.2" + "version": "3.6.8" } }, "nbformat": 4, diff --git a/fast_tsne.R b/fast_tsne.R index 1bbbeba..ff9c939 100644 --- a/fast_tsne.R +++ b/fast_tsne.R @@ -124,7 +124,7 @@ fftRtsne <- function(X, } else if (perplexity==0) { search_k = n_trees*max(perplexity_list)*3 } else { - search_k = n_trees*K*3 + search_k = n_trees*K } } diff --git a/fast_tsne.py b/fast_tsne.py index fbc4c89..563393e 100644 --- a/fast_tsne.py +++ b/fast_tsne.py @@ -18,11 +18,111 @@ def fast_tsne(X, theta=.5, perplexity=30, map_dims=2, max_iter=1000, stop_early_exag_iter=250, K=-1, sigma=-1, nbody_algo='FFT', knn_algo='annoy', mom_switch_iter=250, momentum=.5, final_momentum=.8, learning_rate=200, - early_exag_coeff=12, no_momentum_during_exag=0, n_trees=50, + early_exag_coeff=12, no_momentum_during_exag=False, n_trees=50, search_k=None, start_late_exag_iter=-1, late_exag_coeff=-1, nterms=3, intervals_per_integer=1, min_num_intervals=50, seed=-1, initialization=None, load_affinities=None, - perplexity_list=None, df=1, return_loss=False, nthreads=0): + perplexity_list=None, df=1, return_loss=False, nthreads=None): + """Run t-SNE. This implementation supports exact t-SNE, Barnes-Hut t-SNE and FFT-accelerated + interpolation-based t-SNE (FIt-SNE). This is a Python wrapper to a C++ executable. + + Parameters + ---------- + X: 2D numpy array + Array of observations (n times p) + perplexity: double + Perplexity is used to determine the bandwidth of the Gaussian kernel in the input + space. Default 30. + theta: double + Set to 0 for exact t-SNE. If non-zero, then the code will use either Barnes Hut + or FIt-SNE based on `nbody_algo`. If Barnes Hut, then theta determins the accuracy of + BH approximation. Default 0.5. + map_dims: int + Number of embedding dimensions. Default 2. FIt-SNE supports only 1 or 2 dimensions. + max_iter: int + Number of gradient descent iterations. Default 1000. + nbody_algo: {'Barnes-Hut', 'FFT'} + If theta is nonzero, this determins whether to use FIt-SNE or Barnes Hut approximation. + Default is 'FFT'. + knn_algo: {'vp-tree', 'annoy'} + Use exact nearest neighbours with VP trees (as in BH t-SNE) or approximate nearest neighbors + with Annoy. Default is 'annoy'. + early_exag_coeff: double + Coefficient for early exaggeration. Default 12. + stop_early_exag_iter: int + When to switch off early exaggeration. Default 250. + late_exag_coeff: double + Coefficient for late exaggeration. Set to -1 in order not to use late exaggeration. + Default -1. + start_late_exag_iter: + When to start late exaggeration. Set to -1 in order not to use late exaggeration. + Default -1. + momentum: double + Initial value of momentum. Default 0.5. + final_momentum: double + The value of momentum to use later in the optimisation. Default 0.8. + mom_switch_iter: int + Iteration number to switch from momentum to final_momentum. Default 250. + learning_rate: double + Learning rate. Default 200. + no_mometum_during_exag: boolean + Whether to switch off momentum during the early exaggeration phase (can be useful + for experiments with large exaggeration coefficients). Default is False. + sigma: boolean + The standard deviation of the Gaussian kernel to be used for all points instead of + choosing it adaptively via perplexity. Set to -1 to use perplexity. Default is -1. + K: int + The number of nearest neighbours to use when using fixed sigma instead of perplexity + calibration. Set to -1 when perplexity is used. Default is -1. + nterms: int + If using FIt-SNE, this is the number of interpolation points per sub-interval + intervals_per_integer: double + See min_num_intervals + min_num_intervals: int + The interpolation grid is chosen on each step of the gradient descent. If Y is the current + embedding, let maxloc = ceiling(max(Y.flatten)) and minloc = floor(min(Y.flatten)), i.e. + the points are contained in a [minloc, maxloc]^no_dims box. The number of intervals in each + dimension is either min_num_intervals or ceiling((maxloc-minloc)/intervals_per_integer), + whichever is larger. min_num_intervals must be a positive integer and intervals_per_integer + must be positive real value. Defaults: min_num_intervals=50, intervals_per_integer = 1. + n_trees: int + When using Annoy, the number of search trees to use. Default is 50. + search_k: int + When using Annoy, the number of nodes to inspect during search. Default is 3*perplexity*n_trees + (or K*n_trees when using fixed sigma). + seed: int + Seed for random initialisation. Use -1 to initialise random number generator with current time. + Default -1. + initialization: numpy aray + N x no_dims array to intialize the solution. Default: None. + load_affinities: {'load', 'save', None} + If 'save', input similarities (p_ij) are saved into a file. If 'load', they are loaded from a file + and not recomputed. If None, they are not saved and not loaded. Default is None. + perplexity_list: list + A list of perplexities to used as a perplexity combination. Input affinities are computed with each + perplexity on the list and then averaged. Default is None. + nthreads: int + Number of threads to use. Default is None (i.e. use all available threads). + df: double + Controls the degree of freedom of t-distribution. Must be positive. The actual degree of + freedom is 2*df-1. The standard t-SNE choice of 1 degree of freedom corresponds to df=1. + Large df approximates Gaussian kernel. df<1 corresponds to heavier tails, which can often + resolve substructure in the embedding. See Kobak et al. (2019) for details. Default is 1.0. + return_loss: boolean + If True, the function returns the loss values computed during optimisation + together with the final embedding. If False, only the embedding is returned. + Default is False. + + Returns + ------- + Y: numpy array + The embedding. + loss: numpy array + Loss values computed during optimisation. Only returned if return_loss is True. + """ + + version_number = '1.1.0' + # X should be a numpy array of 64-bit doubles X = np.array(X).astype(float) @@ -39,7 +139,7 @@ def fast_tsne(X, theta=.5, perplexity=30, map_dims=2, max_iter=1000, elif perplexity == 0: search_k = 3 * np.max(perplexity_list) * n_trees else: - search_k = 3 * K * n_trees + search_k = K * n_trees if nbody_algo == 'Barnes-Hut': nbody_algo = 1 @@ -58,6 +158,14 @@ def fast_tsne(X, theta=.5, perplexity=30, map_dims=2, max_iter=1000, else: load_affinities = 0 + if nthreads is None: + nthreads = 0 + + if no_momentum_during_exag: + no_momentum_during_exag = 1 + else: + no_momentum_during_exag = 0 + # write data file with open(os.getcwd() + '/data.dat', 'wb') as f: n, d = X.shape @@ -100,7 +208,7 @@ def fast_tsne(X, theta=.5, perplexity=30, map_dims=2, max_iter=1000, # run t-sne subprocess.call([os.path.dirname(os.path.realpath(__file__)) + - '/bin/fast_tsne', 'data.dat', 'result.dat', '{}'.format(nthreads)]) + '/bin/fast_tsne', version_number, 'data.dat', 'result.dat', '{}'.format(nthreads)]) # read data file with open(os.getcwd()+'/result.dat', 'rb') as f: From 0e9396ee663ec1525e86e46a619918e8b56c7a36 Mon Sep 17 00:00:00 2001 From: Dmitry Kobak Date: Fri, 8 Feb 2019 09:46:55 +0100 Subject: [PATCH 7/9] edited the list of features --- README.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index f1aa370..e4b0b40 100644 --- a/README.md +++ b/README.md @@ -11,11 +11,11 @@ Check out our [preprint](https://arxiv.org/abs/1712.09005) for more details and ## Features Additionally, this implementation includes the following features: * Early exaggeration: In [Linderman and Steinerberger (2017)](https://arxiv.org/abs/1706.02582), we showed that appropriately choosing the early exaggeration coefficient can lead to improved embedding of swissrolls and other synthetic datasets. Early exaggeration is built into all t-SNE implementations; here we highlight its importance as a parameter. -* Late exaggeration: Increasing the exaggeration coefficient late in the optimization process can improve separation of the clusters. For example, [Kobak and Berens (2018)](https://www.biorxiv.org/content/10.1101/453449v1) suggest starting late exaggeration immediately following early exaggeration. -* Perplexity annealing: The perplexity parameter determines the width of the Gaussian kernel, such that small perplexity values uncover local structure while larger values reveal global structure. [Kobak and Berens (2018)](https://www.biorxiv.org/content/10.1101/453449v1) suggest using a range of perplexity values, where the large perplexity is used as initialization for t-SNE with small perplexity, resulting in a multi-scale embedding. -* PCA initialization: As suggested by [Kobak and Berens (2018)](https://www.biorxiv.org/content/10.1101/453449v1), initializing t-SNE with the first two principal components (scaled to have standard deviation 0.0001) results in an embedding which preserves the global structure more effectively than the default random normalization. -* Downsampling-based initialization: +* Late exaggeration: Increasing the exaggeration coefficient late in the optimization process can improve separation of the clusters. [Kobak and Berens (2018)](https://www.biorxiv.org/content/10.1101/453449v1) suggest starting late exaggeration immediately following early exaggeration. +* Initialization: Custom initialization can be provided from Python/Matlab/R. As suggested by [Kobak and Berens (2018)](https://www.biorxiv.org/content/10.1101/453449v1), initializing t-SNE with the first two principal components (scaled to have standard deviation 0.0001) results in an embedding which often preserves the global structure more effectively than the default random normalization. See there for other initialisation tricks. * Variable degrees of freedom: [Kobak et al. (2019)]() show that decreasing the degree of freedom (df) of the t-distribution (resulting in heavier tails) reveals fine structure that is not visible in standard t-SNE. +* Perplexity combination: The perplexity parameter determines the width of the Gaussian kernel, such that small perplexity values uncover local structure while larger values reveal global structure. [Kobak and Berens (2018)](https://www.biorxiv.org/content/10.1101/453449v1) show that using combination of several perplexity values, resulting in a multi-scale embedding, can be useful. +* All optimisation parameters can be controlled from Python/Matlab/R. For example, [Belkina et al. (2018)](https://www.biorxiv.org/content/10.1101/451690v2) highlight the importance of increasing the learning rate when embedding large data sets. ## Installation From b2fdf9c6baf604f2c7df043aa656dd4e4cbb4529 Mon Sep 17 00:00:00 2001 From: linqiaozhi Date: Fri, 8 Feb 2019 10:13:07 -0500 Subject: [PATCH 8/9] Ran notebook with v1.1.0 --- examples/test.ipynb | 1058 +++++++++++++++++++++++++++++++++++++------ 1 file changed, 930 insertions(+), 128 deletions(-) diff --git a/examples/test.ipynb b/examples/test.ipynb index 065f2fb..a969eae 100644 --- a/examples/test.ipynb +++ b/examples/test.ipynb @@ -12,9 +12,7 @@ { "cell_type": "code", "execution_count": 1, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "%matplotlib notebook\n", @@ -31,15 +29,13 @@ { "cell_type": "code", "execution_count": 2, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "/home/dmitry/anaconda3/lib/python3.6/site-packages/h5py/__init__.py:34: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n", + "/Users/george/anaconda3/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n", " from ._conv import register_converters as _register_converters\n", "Using TensorFlow backend.\n" ] @@ -78,9 +74,7 @@ { "cell_type": "code", "execution_count": 3, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "data": { @@ -865,7 +859,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -888,9 +882,7 @@ { "cell_type": "code", "execution_count": 4, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "data": { @@ -974,7 +966,7 @@ " };\n", "\n", " this.imageObj.onunload = function() {\n", - " this.ws.close();\n", + " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", @@ -1469,7 +1461,7 @@ " // Register the callback with on_msg.\n", " comm.on_msg(function(msg) {\n", " //console.log('receiving', msg['content']['data'], msg);\n", - " // Pass the mpl event to the overriden (by mpl) onmessage function.\n", + " // Pass the mpl event to the overridden (by mpl) onmessage function.\n", " ws.onmessage(msg['content']['data'])\n", " });\n", " return ws;\n", @@ -1621,9 +1613,12 @@ " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", - " // select the cell after this one\n", - " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", - " IPython.notebook.select(index + 1);\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", " }\n", "}\n", "\n", @@ -1672,7 +1667,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -1712,10 +1707,8 @@ }, { "cell_type": "code", - "execution_count": 3, - "metadata": { - "collapsed": true - }, + "execution_count": 5, + "metadata": {}, "outputs": [], "source": [ "# Subsampling \n", @@ -1735,9 +1728,7 @@ { "cell_type": "code", "execution_count": 6, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "data": { @@ -1821,7 +1812,7 @@ " };\n", "\n", " this.imageObj.onunload = function() {\n", - " this.ws.close();\n", + " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", @@ -2316,7 +2307,7 @@ " // Register the callback with on_msg.\n", " comm.on_msg(function(msg) {\n", " //console.log('receiving', msg['content']['data'], msg);\n", - " // Pass the mpl event to the overriden (by mpl) onmessage function.\n", + " // Pass the mpl event to the overridden (by mpl) onmessage function.\n", " ws.onmessage(msg['content']['data'])\n", " });\n", " return ws;\n", @@ -2468,9 +2459,12 @@ " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", - " // select the cell after this one\n", - " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", - " IPython.notebook.select(index + 1);\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", " }\n", "}\n", "\n", @@ -2519,7 +2513,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -2557,9 +2551,7 @@ { "cell_type": "code", "execution_count": 7, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "data": { @@ -2643,7 +2635,7 @@ " };\n", "\n", " this.imageObj.onunload = function() {\n", - " this.ws.close();\n", + " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", @@ -3138,7 +3130,7 @@ " // Register the callback with on_msg.\n", " comm.on_msg(function(msg) {\n", " //console.log('receiving', msg['content']['data'], msg);\n", - " // Pass the mpl event to the overriden (by mpl) onmessage function.\n", + " // Pass the mpl event to the overridden (by mpl) onmessage function.\n", " ws.onmessage(msg['content']['data'])\n", " });\n", " return ws;\n", @@ -3290,9 +3282,12 @@ " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", - " // select the cell after this one\n", - " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", - " IPython.notebook.select(index + 1);\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", " }\n", "}\n", "\n", @@ -3341,7 +3336,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -3378,10 +3373,8 @@ }, { "cell_type": "code", - "execution_count": 6, - "metadata": { - "collapsed": false - }, + "execution_count": 8, + "metadata": {}, "outputs": [ { "data": { @@ -3465,7 +3458,7 @@ " };\n", "\n", " this.imageObj.onunload = function() {\n", - " this.ws.close();\n", + " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", @@ -3960,7 +3953,7 @@ " // Register the callback with on_msg.\n", " comm.on_msg(function(msg) {\n", " //console.log('receiving', msg['content']['data'], msg);\n", - " // Pass the mpl event to the overriden (by mpl) onmessage function.\n", + " // Pass the mpl event to the overridden (by mpl) onmessage function.\n", " ws.onmessage(msg['content']['data'])\n", " });\n", " return ws;\n", @@ -4112,9 +4105,12 @@ " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", - " // select the cell after this one\n", - " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", - " IPython.notebook.select(index + 1);\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", " }\n", "}\n", "\n", @@ -4163,7 +4159,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -4200,10 +4196,8 @@ }, { "cell_type": "code", - "execution_count": 11, - "metadata": { - "collapsed": false - }, + "execution_count": 9, + "metadata": {}, "outputs": [ { "data": { @@ -4287,7 +4281,7 @@ " };\n", "\n", " this.imageObj.onunload = function() {\n", - " this.ws.close();\n", + " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", @@ -4782,7 +4776,7 @@ " // Register the callback with on_msg.\n", " comm.on_msg(function(msg) {\n", " //console.log('receiving', msg['content']['data'], msg);\n", - " // Pass the mpl event to the overriden (by mpl) onmessage function.\n", + " // Pass the mpl event to the overridden (by mpl) onmessage function.\n", " ws.onmessage(msg['content']['data'])\n", " });\n", " return ws;\n", @@ -4934,9 +4928,12 @@ " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", - " // select the cell after this one\n", - " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", - " IPython.notebook.select(index + 1);\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", " }\n", "}\n", "\n", @@ -4985,7 +4982,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -5021,10 +5018,8 @@ }, { "cell_type": "code", - "execution_count": 9, - "metadata": { - "collapsed": false - }, + "execution_count": 10, + "metadata": {}, "outputs": [ { "data": { @@ -5108,7 +5103,7 @@ " };\n", "\n", " this.imageObj.onunload = function() {\n", - " this.ws.close();\n", + " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", @@ -5603,7 +5598,7 @@ " // Register the callback with on_msg.\n", " comm.on_msg(function(msg) {\n", " //console.log('receiving', msg['content']['data'], msg);\n", - " // Pass the mpl event to the overriden (by mpl) onmessage function.\n", + " // Pass the mpl event to the overridden (by mpl) onmessage function.\n", " ws.onmessage(msg['content']['data'])\n", " });\n", " return ws;\n", @@ -5755,9 +5750,12 @@ " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", - " // select the cell after this one\n", - " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", - " IPython.notebook.select(index + 1);\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", " }\n", "}\n", "\n", @@ -5806,7 +5804,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -5841,10 +5839,8 @@ }, { "cell_type": "code", - "execution_count": 12, - "metadata": { - "collapsed": false - }, + "execution_count": 11, + "metadata": {}, "outputs": [ { "data": { @@ -5928,7 +5924,7 @@ " };\n", "\n", " this.imageObj.onunload = function() {\n", - " this.ws.close();\n", + " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", @@ -6423,7 +6419,7 @@ " // Register the callback with on_msg.\n", " comm.on_msg(function(msg) {\n", " //console.log('receiving', msg['content']['data'], msg);\n", - " // Pass the mpl event to the overriden (by mpl) onmessage function.\n", + " // Pass the mpl event to the overridden (by mpl) onmessage function.\n", " ws.onmessage(msg['content']['data'])\n", " });\n", " return ws;\n", @@ -6575,9 +6571,12 @@ " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", - " // select the cell after this one\n", - " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", - " IPython.notebook.select(index + 1);\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", " }\n", "}\n", "\n", @@ -6626,7 +6625,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -6717,7 +6716,7 @@ " };\n", "\n", " this.imageObj.onunload = function() {\n", - " this.ws.close();\n", + " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", @@ -7212,7 +7211,7 @@ " // Register the callback with on_msg.\n", " comm.on_msg(function(msg) {\n", " //console.log('receiving', msg['content']['data'], msg);\n", - " // Pass the mpl event to the overriden (by mpl) onmessage function.\n", + " // Pass the mpl event to the overridden (by mpl) onmessage function.\n", " ws.onmessage(msg['content']['data'])\n", " });\n", " return ws;\n", @@ -7364,9 +7363,12 @@ " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", - " // select the cell after this one\n", - " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", - " IPython.notebook.select(index + 1);\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", " }\n", "}\n", "\n", @@ -7415,7 +7417,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -7472,10 +7474,8 @@ }, { "cell_type": "code", - "execution_count": 21, - "metadata": { - "collapsed": false - }, + "execution_count": 12, + "metadata": {}, "outputs": [ { "data": { @@ -7559,7 +7559,7 @@ " };\n", "\n", " this.imageObj.onunload = function() {\n", - " this.ws.close();\n", + " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", @@ -8054,7 +8054,7 @@ " // Register the callback with on_msg.\n", " comm.on_msg(function(msg) {\n", " //console.log('receiving', msg['content']['data'], msg);\n", - " // Pass the mpl event to the overriden (by mpl) onmessage function.\n", + " // Pass the mpl event to the overridden (by mpl) onmessage function.\n", " ws.onmessage(msg['content']['data'])\n", " });\n", " return ws;\n", @@ -8206,9 +8206,12 @@ " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", - " // select the cell after this one\n", - " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", - " IPython.notebook.select(index + 1);\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", " }\n", "}\n", "\n", @@ -8257,7 +8260,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -8293,10 +8296,8 @@ }, { "cell_type": "code", - "execution_count": 10, - "metadata": { - "collapsed": false - }, + "execution_count": 13, + "metadata": {}, "outputs": [ { "data": { @@ -8380,7 +8381,7 @@ " };\n", "\n", " this.imageObj.onunload = function() {\n", - " this.ws.close();\n", + " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", @@ -8875,7 +8876,7 @@ " // Register the callback with on_msg.\n", " comm.on_msg(function(msg) {\n", " //console.log('receiving', msg['content']['data'], msg);\n", - " // Pass the mpl event to the overriden (by mpl) onmessage function.\n", + " // Pass the mpl event to the overridden (by mpl) onmessage function.\n", " ws.onmessage(msg['content']['data'])\n", " });\n", " return ws;\n", @@ -9027,9 +9028,12 @@ " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", - " // select the cell after this one\n", - " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", - " IPython.notebook.select(index + 1);\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", " }\n", "}\n", "\n", @@ -9078,7 +9082,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -9114,25 +9118,823 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ - "Z1 = fast_tsne(X50, perplexity=50, learning_rate=1000, initialization=PCAinit)\n", - "\n", - "Z2 = fast_tsne(X50, perplexity=50, learning_rate=1000, initialization=PCAinit, df=.5)\n", - "\n", - "plt.figure(figsize=(8,4))\n", - "plt.subplot(121)\n", - "plt.scatter(Z1[:,0], Z1[:,1], c=col[y], s=1)\n", - "plt.title('1/(1+x^2) kernel')\n", - "plt.subplot(122)\n", - "plt.scatter(Z2[:,0], Z2[:,1], c=col[y], s=1)\n", - "plt.title('1/(1+x^2/0.5)^0.5 kernel')\n", - "plt.tight_layout()" - ] + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support.' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " if (mpl.ratio != 1) {\n", + " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", + " }\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " fig.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var backingStore = this.context.backingStorePixelRatio ||\n", + "\tthis.context.webkitBackingStorePixelRatio ||\n", + "\tthis.context.mozBackingStorePixelRatio ||\n", + "\tthis.context.msBackingStorePixelRatio ||\n", + "\tthis.context.oBackingStorePixelRatio ||\n", + "\tthis.context.backingStorePixelRatio || 1;\n", + "\n", + " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width * mpl.ratio);\n", + " canvas.attr('height', height * mpl.ratio);\n", + " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
')\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('