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 <dmitry.kobak@neuro.fchampalimaud.org>

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 {

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

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;
         }