LCOV - code coverage report
Current view: top level - external/barnes_hut_sne - tsne.hpp (source / functions) Hit Total Coverage
Test: clean.info Lines: 188 271 69.4 %
Date: 2013-05-24 Functions: 6 9 66.7 %
Branches: 175 328 53.4 %

           Branch data     Line data    Source code
       1                 :            : /**
       2                 :            :  * Copyright (c) 2013, Laurens van der Maaten (Delft University of Technology)
       3                 :            :  * All rights reserved.
       4                 :            :  *
       5                 :            :  * Redistribution and use in source and binary forms, with or without
       6                 :            :  * modification, are permitted provided that the following conditions are met:
       7                 :            :  * 1. Redistributions of source code must retain the above copyright
       8                 :            :  *    notice, this list of conditions and the following disclaimer.
       9                 :            :  * 2. Redistributions in binary form must reproduce the above copyright
      10                 :            :  *    notice, this list of conditions and the following disclaimer in the
      11                 :            :  *    documentation and/or other materials provided with the distribution.
      12                 :            :  * 3. All advertising materials mentioning features or use of this software
      13                 :            :  *    must display the following acknowledgement:
      14                 :            :  *    This product includes software developed by the Delft University of Technology.
      15                 :            :  * 4. Neither the name of the Delft University of Technology nor the names of
      16                 :            :  *    its contributors may be used to endorse or promote products derived from
      17                 :            :  *    this software without specific prior written permission.
      18                 :            :  *
      19                 :            :  * THIS SOFTWARE IS PROVIDED BY LAURENS VAN DER MAATEN ''AS IS'' AND ANY EXPRESS
      20                 :            :  * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
      21                 :            :  * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
      22                 :            :  * EVENT SHALL LAURENS VAN DER MAATEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
      23                 :            :  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
      24                 :            :  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
      25                 :            :  * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
      26                 :            :  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
      27                 :            :  * IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
      28                 :            :  * OF SUCH DAMAGE.
      29                 :            :  *
      30                 :            :  */
      31                 :            : 
      32                 :            : #ifndef TSNE_H
      33                 :            : #define TSNE_H
      34                 :            : 
      35                 :            : /* Tapkee includes */
      36                 :            : #include <tapkee/utils/logging.hpp>
      37                 :            : #include <tapkee/utils/time.hpp>
      38                 :            : #include <tapkee/external/barnes_hut_sne/quadtree.hpp>
      39                 :            : #include <tapkee/external/barnes_hut_sne/vptree.hpp>
      40                 :            : /* End of Tapkee includes */
      41                 :            : 
      42                 :            : #include <math.h>
      43                 :            : #include <float.h>
      44                 :            : #include <stdlib.h>
      45                 :            : #include <stdio.h>
      46                 :            : #include <cstring>
      47                 :            : #include <time.h>
      48                 :            : 
      49                 :            : //! Namespace containing implementation of t-SNE algorithm
      50                 :            : namespace tsne
      51                 :            : {
      52                 :            : 
      53 [ +  + ][ +  + ]:     200000 : static inline double sign(double x) { return (x == .0 ? .0 : (x < .0 ? -1.0 : 1.0)); }
      54                 :            : 
      55                 :            : class TSNE
      56                 :            : {    
      57                 :            : public:
      58                 :          1 :         void run(double* X, int N, int D, double* Y, int no_dims, double perplexity, double theta)
      59                 :            :         {
      60                 :            :                 // Determine whether we are using an exact algorithm
      61                 :          1 :                 bool exact = (theta == .0) ? true : false;
      62         [ -  + ]:          1 :                 if (exact) 
      63 [ #  # ][ #  # ]:          0 :                         tapkee::LoggingSingleton::instance().message_info("Using exact t-SNE algorithm");
         [ #  # ][ #  # ]
      64                 :            :                 else
      65 [ +  - ][ +  - ]:          1 :                         tapkee::LoggingSingleton::instance().message_info("Using Barnes-Hut-SNE algorithm");
         [ +  - ][ +  - ]
      66                 :            :                 
      67                 :            :                 // Set learning parameters
      68                 :          1 :                 float total_time = .0;
      69                 :            :                 clock_t start, end;
      70                 :          1 :                 int max_iter = 1000, stop_lying_iter = 250, mom_switch_iter = 250;
      71                 :          1 :                 double momentum = .5, final_momentum = .8;
      72                 :          1 :                 double eta = 200.0;
      73                 :            :                 
      74                 :            :                 // Allocate some memory
      75                 :          1 :                 double* dY    = (double*) malloc(N * no_dims * sizeof(double));
      76                 :          1 :                 double* uY    = (double*) malloc(N * no_dims * sizeof(double));
      77                 :          1 :                 double* gains = (double*) malloc(N * no_dims * sizeof(double));
      78 [ +  - ][ +  - ]:          1 :                 if(dY == NULL || uY == NULL || gains == NULL) { printf("Memory allocation failed!\n"); exit(1); }
                 [ -  + ]
      79         [ +  + ]:        101 :                 for(int i = 0; i < N * no_dims; i++)    uY[i] =  .0;
      80         [ +  + ]:        101 :                 for(int i = 0; i < N * no_dims; i++) gains[i] = 1.0;
      81                 :            :                 
      82                 :            :                 // Normalize input data (to prevent numerical problems)
      83                 :          1 :                 double* P=NULL; int* row_P; int* col_P; double* val_P;
      84                 :            :                 {
      85 [ +  - ][ +  - ]:          1 :                         tapkee::tapkee_internal::timed_context context("Input similarities computation");
                 [ +  - ]
      86                 :          1 :                         start = clock();
      87         [ +  - ]:          1 :                         zeroMean(X, N, D);
      88                 :          1 :                         double max_X = .0;
      89         [ +  + ]:        151 :                         for(int i = 0; i < N * D; i++) {
      90         [ +  + ]:        150 :                                 if(X[i] > max_X) max_X = X[i];
      91                 :            :                         }
      92         [ +  + ]:        151 :                         for(int i = 0; i < N * D; i++) X[i] /= max_X;
      93                 :            :                         
      94                 :            :                         // Compute input similarities for exact t-SNE
      95         [ -  + ]:          1 :                         if(exact) {
      96                 :            :                                 
      97                 :            :                                 // Compute similarities
      98                 :          0 :                                 P = (double*) malloc(N * N * sizeof(double));
      99 [ #  # ][ #  # ]:          0 :                                 if(P == NULL) { printf("Memory allocation failed!\n"); exit(1); }
     100         [ #  # ]:          0 :                                 computeGaussianPerplexity(X, N, D, P, perplexity);
     101                 :            :                         
     102                 :            :                                 // Symmetrize input similarities
     103         [ #  # ]:          0 :                                 for(int n = 0; n < N; n++) {
     104         [ #  # ]:          0 :                                         for(int m = n + 1; m < N; m++) {
     105                 :          0 :                                                 P[n * N + m] += P[m * N + n];
     106                 :          0 :                                                 P[m * N + n]  = P[n * N + m];
     107                 :            :                                         }
     108                 :            :                                 }
     109                 :          0 :                                 double sum_P = .0;
     110         [ #  # ]:          0 :                                 for(int i = 0; i < N * N; i++) sum_P += P[i];
     111         [ #  # ]:          0 :                                 for(int i = 0; i < N * N; i++) P[i] /= sum_P;
     112                 :            :                         }
     113                 :            :                         
     114                 :            :                         // Compute input similarities for approximate t-SNE
     115                 :            :                         else {
     116                 :            :                         
     117                 :            :                                 // Compute asymmetric pairwise input similarities
     118         [ +  - ]:          1 :                                 computeGaussianPerplexity(X, N, D, &row_P, &col_P, &val_P, perplexity, (int) (3 * perplexity));
     119                 :            :                                 
     120                 :            :                                 // Symmetrize input similarities
     121         [ +  - ]:          1 :                                 symmetrizeMatrix(&row_P, &col_P, &val_P, N);
     122                 :          1 :                                 double sum_P = .0;
     123         [ +  + ]:       1745 :                                 for(int i = 0; i < row_P[N]; i++) sum_P += val_P[i];
     124         [ +  + ]:       1745 :                                 for(int i = 0; i < row_P[N]; i++) val_P[i] /= sum_P;
     125                 :            :                         }
     126                 :            :                         
     127                 :            :                         // Lie about the P-values
     128 [ -  + ][ #  # ]:          1 :                         if(exact) { for(int i = 0; i < N * N; i++)        P[i] *= 12.0; }
     129         [ +  + ]:       1745 :                         else {      for(int i = 0; i < row_P[N]; i++) val_P[i] *= 12.0; }
     130                 :            : 
     131                 :            :                         // Initialize solution (randomly)
     132         [ +  + ]:        101 :                         for(int i = 0; i < N * no_dims; i++) Y[i] = tapkee::gaussian_random() * .0001;
     133                 :            :                 }
     134                 :            : 
     135                 :            :                 {
     136 [ +  - ][ +  - ]:          1 :                         tapkee::tapkee_internal::timed_context context("Main t-SNE loop");
                 [ +  - ]
     137         [ +  + ]:       1001 :                         for(int iter = 0; iter < max_iter; iter++) {
     138                 :            :                                 
     139                 :            :                                 // Compute (approximate) gradient
     140 [ -  + ][ #  # ]:       1000 :                                 if(exact) computeExactGradient(P, Y, N, no_dims, dY);
     141         [ +  - ]:       1000 :                                 else computeGradient(P, row_P, col_P, val_P, Y, N, no_dims, dY, theta);
     142                 :            :                                 
     143                 :            :                                 // Update gains
     144 [ +  + ][ +  + ]:     101000 :                                 for(int i = 0; i < N * no_dims; i++) gains[i] = (sign(dY[i]) != sign(uY[i])) ? (gains[i] + .2) : (gains[i] * .8);
     145 [ +  + ][ +  + ]:     101000 :                                 for(int i = 0; i < N * no_dims; i++) if(gains[i] < .01) gains[i] = .01;
     146                 :            :                                         
     147                 :            :                                 // Perform gradient update (with momentum and gains)
     148         [ +  + ]:     101000 :                                 for(int i = 0; i < N * no_dims; i++) uY[i] = momentum * uY[i] - eta * gains[i] * dY[i];
     149         [ +  + ]:     101000 :                                 for(int i = 0; i < N * no_dims; i++)  Y[i] = Y[i] + uY[i];
     150                 :            :                                 
     151                 :            :                                 // Make solution zero-mean
     152         [ +  - ]:       1000 :                                 zeroMean(Y, N, no_dims);
     153                 :            :                                 
     154                 :            :                                 // Stop lying about the P-values after a while, and switch momentum
     155         [ +  + ]:       1000 :                                 if(iter == stop_lying_iter) {
     156 [ -  + ][ #  # ]:          1 :                                         if(exact) { for(int i = 0; i < N * N; i++)        P[i] /= 12.0; }
     157         [ +  + ]:       1745 :                                         else      { for(int i = 0; i < row_P[N]; i++) val_P[i] /= 12.0; }
     158                 :            :                                 }
     159         [ +  + ]:       1000 :                                 if(iter == mom_switch_iter) momentum = final_momentum;
     160                 :            :                                 
     161                 :            :                                 // Print out progress
     162                 :            :                                 /*
     163                 :            :                                 if((iter > 0) && ((iter % 50 == 0) || (iter == max_iter - 1))) {
     164                 :            :                                         end = clock();
     165                 :            :                                         double C = .0;
     166                 :            :                                         if(exact) C = evaluateError(P, Y, N);
     167                 :            :                                         else      C = evaluateError(row_P, col_P, val_P, Y, N, theta);  // doing approximate computation here!
     168                 :            :                                         if(iter == 0) 
     169                 :            :                                         {
     170                 :            :                                                 //printf("Iteration %d: error is %f\n", iter + 1, C);
     171                 :            :                                         }
     172                 :            :                                         else 
     173                 :            :                                         {
     174                 :            :                                                 total_time += (float) (end - start) / CLOCKS_PER_SEC;
     175                 :            :                                                 //printf("Iteration %d: error is %f (50 iterations in %4.2f seconds)\n", iter, C, (float) (end - start) / CLOCKS_PER_SEC);
     176                 :            :                                         }
     177                 :            :                                         start = clock();
     178                 :            :                                 }
     179                 :            :                                 */
     180                 :            :                         }
     181                 :          1 :                         end = clock(); total_time += (float) (end - start) / CLOCKS_PER_SEC;
     182                 :            :                         
     183                 :            :                         // Clean up memory
     184                 :          1 :                         free(dY);
     185                 :          1 :                         free(uY);
     186                 :          1 :                         free(gains);
     187         [ -  + ]:          1 :                         if(exact) free(P);
     188                 :            :                         else {
     189                 :          1 :                                 free(row_P); row_P = NULL;
     190                 :          1 :                                 free(col_P); col_P = NULL;
     191                 :          1 :                                 free(val_P); val_P = NULL;
     192                 :          1 :                         }
     193                 :            :                 }
     194                 :          1 :         }
     195                 :            : 
     196                 :          1 :         void symmetrizeMatrix(int** _row_P, int** _col_P, double** _val_P, int N)
     197                 :            :         { 
     198                 :            :                 // Get sparse matrix
     199                 :          1 :                 int* row_P = *_row_P;
     200                 :          1 :                 int* col_P = *_col_P;
     201                 :          1 :                 double* val_P = *_val_P;
     202                 :            : 
     203                 :            :                 // Count number of elements and row counts of symmetric matrix
     204                 :          1 :                 int* row_counts = (int*) calloc(N, sizeof(int));
     205         [ -  + ]:          1 :                 if(row_counts == NULL) { printf("Memory allocation failed!\n"); exit(1); }
     206         [ +  + ]:         51 :                 for(int n = 0; n < N; n++) {
     207         [ +  + ]:       1550 :                         for(int i = row_P[n]; i < row_P[n + 1]; i++) {
     208                 :            :                                 // Check whether element (col_P[i], n) is present
     209                 :       1500 :                                 bool present = false;
     210         [ +  + ]:      46500 :                                 for(int m = row_P[col_P[i]]; m < row_P[col_P[i] + 1]; m++) {
     211         [ +  + ]:      45000 :                                         if(col_P[m] == n) present = true;
     212                 :            :                                 }
     213         [ +  + ]:       1500 :                                 if(present) row_counts[n]++;
     214                 :            :                                 else {
     215                 :        244 :                                         row_counts[n]++;
     216                 :        244 :                                         row_counts[col_P[i]]++;
     217                 :            :                                 }
     218                 :            :                         }
     219                 :            :                 }
     220                 :          1 :                 int no_elem = 0;
     221         [ +  + ]:         51 :                 for(int n = 0; n < N; n++) no_elem += row_counts[n];
     222                 :            :                 
     223                 :            :                 // Allocate memory for symmetrized matrix
     224                 :          1 :                 int*    sym_row_P = (int*)    malloc((N + 1) * sizeof(int));
     225                 :          1 :                 int*    sym_col_P = (int*)    malloc(no_elem * sizeof(int));
     226                 :          1 :                 double* sym_val_P = (double*) malloc(no_elem * sizeof(double));
     227 [ +  - ][ +  - ]:          1 :                 if(sym_row_P == NULL || sym_col_P == NULL || sym_val_P == NULL) { printf("Memory allocation failed!\n"); exit(1); }
                 [ -  + ]
     228                 :            :                 
     229                 :            :                 // Construct new row indices for symmetric matrix
     230                 :          1 :                 sym_row_P[0] = 0;
     231         [ +  + ]:         51 :                 for(int n = 0; n < N; n++) sym_row_P[n + 1] = sym_row_P[n] + row_counts[n];
     232                 :            :                 
     233                 :            :                 // Fill the result matrix
     234                 :          1 :                 int* offset = (int*) calloc(N, sizeof(int));
     235         [ -  + ]:          1 :                 if(offset == NULL) { printf("Memory allocation failed!\n"); exit(1); }
     236         [ +  + ]:         51 :                 for(int n = 0; n < N; n++) {
     237         [ +  + ]:       1550 :                         for(int i = row_P[n]; i < row_P[n + 1]; i++) {                                  // considering element(n, col_P[i])
     238                 :            :                                 
     239                 :            :                                 // Check whether element (col_P[i], n) is present
     240                 :       1500 :                                 bool present = false;
     241         [ +  + ]:      46500 :                                 for(int m = row_P[col_P[i]]; m < row_P[col_P[i] + 1]; m++) {
     242         [ +  + ]:      45000 :                                         if(col_P[m] == n) {
     243                 :       1256 :                                                 present = true;
     244         [ +  + ]:       1256 :                                                 if(n <= col_P[i]) {                                                 // make sure we do not add elements twice
     245                 :        628 :                                                         sym_col_P[sym_row_P[n]        + offset[n]]        = col_P[i];
     246                 :        628 :                                                         sym_col_P[sym_row_P[col_P[i]] + offset[col_P[i]]] = n;
     247                 :        628 :                                                         sym_val_P[sym_row_P[n]        + offset[n]]        = val_P[i] + val_P[m];
     248                 :        628 :                                                         sym_val_P[sym_row_P[col_P[i]] + offset[col_P[i]]] = val_P[i] + val_P[m];
     249                 :            :                                                 }
     250                 :            :                                         }
     251                 :            :                                 }
     252                 :            :                                 
     253                 :            :                                 // If (col_P[i], n) is not present, there is no addition involved
     254         [ +  + ]:       1500 :                                 if(!present) {
     255                 :        244 :                                         sym_col_P[sym_row_P[n]        + offset[n]]        = col_P[i];
     256                 :        244 :                                         sym_col_P[sym_row_P[col_P[i]] + offset[col_P[i]]] = n;
     257                 :        244 :                                         sym_val_P[sym_row_P[n]        + offset[n]]        = val_P[i];
     258                 :        244 :                                         sym_val_P[sym_row_P[col_P[i]] + offset[col_P[i]]] = val_P[i];
     259                 :            :                                 }
     260                 :            :                                 
     261                 :            :                                 // Update offsets
     262 [ +  + ][ +  - ]:       1500 :                                 if(!present || (present && n <= col_P[i])) {
                 [ +  + ]
     263                 :        872 :                                         offset[n]++;
     264         [ +  - ]:        872 :                                         if(col_P[i] != n) offset[col_P[i]]++;               
     265                 :            :                                 }
     266                 :            :                         }
     267                 :            :                 }
     268                 :            :                 
     269                 :            :                 // Divide the result by two
     270         [ +  + ]:       1745 :                 for(int i = 0; i < no_elem; i++) sym_val_P[i] /= 2.0;
     271                 :            :                 
     272                 :            :                 // Return symmetrized matrices
     273                 :          1 :                 free(*_row_P); *_row_P = sym_row_P;
     274                 :          1 :                 free(*_col_P); *_col_P = sym_col_P;
     275                 :          1 :                 free(*_val_P); *_val_P = sym_val_P;
     276                 :            :                 
     277                 :            :                 // Free up some memery
     278                 :          1 :                 free(offset); offset = NULL;
     279                 :          1 :                 free(row_counts); row_counts  = NULL;
     280                 :          1 :         }
     281                 :            :     
     282                 :            : private:
     283                 :            : 
     284                 :       1000 :         void computeGradient(double* /*P*/, int* inp_row_P, int* inp_col_P, double* inp_val_P, double* Y, int N, int D, double* dC, double theta)
     285                 :            :         {
     286                 :            :                 // Construct quadtree on current map
     287         [ +  - ]:       1000 :                 QuadTree* tree = new QuadTree(Y, N);
     288                 :            :                 
     289                 :            :                 // Compute all terms required for t-SNE gradient
     290                 :       1000 :                 double sum_Q = .0;
     291                 :       1000 :                 double* pos_f = (double*) calloc(N * D, sizeof(double));
     292                 :       1000 :                 double* neg_f = (double*) calloc(N * D, sizeof(double));
     293 [ +  - ][ -  + ]:       1000 :                 if(pos_f == NULL || neg_f == NULL) { printf("Memory allocation failed!\n"); exit(1); }
     294                 :       1000 :                 tree->computeEdgeForces(inp_row_P, inp_col_P, inp_val_P, N, pos_f);
     295         [ +  + ]:      51000 :                 for(int n = 0; n < N; n++) tree->computeNonEdgeForces(n, theta, neg_f + n * D, &sum_Q);
     296                 :            :                 
     297                 :            :                 // Compute final t-SNE gradient
     298         [ +  + ]:     101000 :                 for(int i = 0; i < N * D; i++) {
     299                 :     100000 :                         dC[i] = pos_f[i] - (neg_f[i] / sum_Q);
     300                 :            :                 }
     301                 :       1000 :                 free(pos_f);
     302                 :       1000 :                 free(neg_f);
     303         [ +  - ]:       1000 :                 delete tree;
     304                 :       1000 :         }
     305                 :            : 
     306                 :          0 :         void computeExactGradient(double* P, double* Y, int N, int D, double* dC)
     307                 :            :         {
     308                 :            :                 // Make sure the current gradient contains zeros
     309         [ #  # ]:          0 :                 for(int i = 0; i < N * D; i++) dC[i] = 0.0;
     310                 :            :                 
     311                 :            :                 // Compute the squared Euclidean distance matrix
     312                 :          0 :                 double* DD = (double*) malloc(N * N * sizeof(double));
     313         [ #  # ]:          0 :                 if(DD == NULL) { printf("Memory allocation failed!\n"); exit(1); }
     314                 :          0 :                 computeSquaredEuclideanDistance(Y, N, D, DD);
     315                 :            :                 
     316                 :            :                 // Compute Q-matrix and normalization sum
     317                 :          0 :                 double* Q    = (double*) malloc(N * N * sizeof(double));
     318         [ #  # ]:          0 :                 if(Q == NULL) { printf("Memory allocation failed!\n"); exit(1); }
     319                 :          0 :                 double sum_Q = .0;
     320         [ #  # ]:          0 :                 for(int n = 0; n < N; n++) {
     321         [ #  # ]:          0 :                         for(int m = 0; m < N; m++) {
     322         [ #  # ]:          0 :                                 if(n != m) {
     323                 :          0 :                                         Q[n * N + m] = 1 / (1 + DD[n * N + m]);
     324                 :          0 :                                         sum_Q += Q[n * N + m];
     325                 :            :                                 }
     326                 :            :                         }
     327                 :            :                 }
     328                 :            :                 
     329                 :            :                 // Perform the computation of the gradient
     330         [ #  # ]:          0 :                 for(int n = 0; n < N; n++) {
     331         [ #  # ]:          0 :                         for(int m = 0; m < N; m++) {
     332         [ #  # ]:          0 :                                 if(n != m) {
     333                 :          0 :                                         double mult = (P[n * N + m] - (Q[n * N + m] / sum_Q)) * Q[n * N + m];
     334         [ #  # ]:          0 :                                         for(int d = 0; d < D; d++) {
     335                 :          0 :                                                 dC[n * D + d] += (Y[n * D + d] - Y[m * D + d]) * mult;
     336                 :            :                                         }
     337                 :            :                                 }
     338                 :            :                         }
     339                 :            :                 }
     340                 :            :                 
     341                 :            :                 // Free memory
     342                 :          0 :                 free(DD); DD = NULL;
     343                 :          0 :                 free(Q);  Q  = NULL;
     344                 :          0 :         }
     345                 :            : 
     346                 :            :         double evaluateError(double* P, double* Y, int N)
     347                 :            :         { 
     348                 :            :                 // Compute the squared Euclidean distance matrix
     349                 :            :                 double* DD = (double*) malloc(N * N * sizeof(double));
     350                 :            :                 double* Q = (double*) malloc(N * N * sizeof(double));
     351                 :            :                 if(DD == NULL || Q == NULL) { printf("Memory allocation failed!\n"); exit(1); }
     352                 :            :                 computeSquaredEuclideanDistance(Y, N, 2, DD);
     353                 :            :                 
     354                 :            :                 // Compute Q-matrix and normalization sum
     355                 :            :                 double sum_Q = DBL_MIN;
     356                 :            :                 for(int n = 0; n < N; n++) {
     357                 :            :                         for(int m = 0; m < N; m++) {
     358                 :            :                                 if(n != m) {
     359                 :            :                                         Q[n * N + m] = 1 / (1 + DD[n * N + m]);
     360                 :            :                                         sum_Q += Q[n * N + m];
     361                 :            :                                 }
     362                 :            :                                 else Q[n * N + m] = DBL_MIN;
     363                 :            :                         }
     364                 :            :                 }
     365                 :            :                 for(int i = 0; i < N * N; i++) Q[i] /= sum_Q;
     366                 :            :                 
     367                 :            :                 // Sum t-SNE error
     368                 :            :                 double C = .0;
     369                 :            :                 for(int n = 0; n < N; n++) {
     370                 :            :                         for(int m = 0; m < N; m++) {
     371                 :            :                                 C += P[n * N + m] * log((P[n * N + m] + 1e-9) / (Q[n * N + m] + 1e-9));
     372                 :            :                         }
     373                 :            :                 }
     374                 :            :                 
     375                 :            :                 // Clean up memory
     376                 :            :                 free(DD);
     377                 :            :                 free(Q);
     378                 :            :                 return C;
     379                 :            :         }
     380                 :            : 
     381                 :            :         double evaluateError(int* row_P, int* col_P, double* val_P, double* Y, int N, double theta)
     382                 :            :         {
     383                 :            :                 // Get estimate of normalization term
     384                 :            :                 const int QT_NO_DIMS = 2;
     385                 :            :                 QuadTree* tree = new QuadTree(Y, N);
     386                 :            :                 double buff[QT_NO_DIMS] = {.0, .0};
     387                 :            :                 double sum_Q = .0;
     388                 :            :                 for(int n = 0; n < N; n++) tree->computeNonEdgeForces(n, theta, buff, &sum_Q);
     389                 :            :                 
     390                 :            :                 // Loop over all edges to compute t-SNE error
     391                 :            :                 int ind1, ind2;
     392                 :            :                 double C = .0, Q;
     393                 :            :                 for(int n = 0; n < N; n++) {
     394                 :            :                         ind1 = n * QT_NO_DIMS;
     395                 :            :                         for(int i = row_P[n]; i < row_P[n + 1]; i++) {
     396                 :            :                                 Q = .0;
     397                 :            :                                 ind2 = col_P[i] * QT_NO_DIMS;
     398                 :            :                                 for(int d = 0; d < QT_NO_DIMS; d++) buff[d]  = Y[ind1 + d];
     399                 :            :                                 for(int d = 0; d < QT_NO_DIMS; d++) buff[d] -= Y[ind2 + d];
     400                 :            :                                 for(int d = 0; d < QT_NO_DIMS; d++) Q += buff[d] * buff[d];
     401                 :            :                                 Q = (1.0 / (1.0 + Q)) / sum_Q;
     402                 :            :                                 C += val_P[i] * log((val_P[i] + FLT_MIN) / (Q + FLT_MIN));
     403                 :            :                         }
     404                 :            :                 }
     405                 :            :                 return C;
     406                 :            :         }
     407                 :            : 
     408                 :       1001 :         void zeroMean(double* X, int N, int D)
     409                 :            :         {
     410                 :            :                 // Compute data mean
     411                 :       1001 :                 double* mean = (double*) calloc(D, sizeof(double));
     412         [ -  + ]:       1001 :                 if(mean == NULL) { printf("Memory allocation failed!\n"); exit(1); }
     413         [ +  + ]:      51051 :                 for(int n = 0; n < N; n++) {
     414         [ +  + ]:     150200 :                         for(int d = 0; d < D; d++) {
     415                 :     100150 :                                 mean[d] += X[n * D + d];
     416                 :            :                         }
     417                 :            :                 }
     418         [ +  + ]:       3004 :                 for(int d = 0; d < D; d++) {
     419                 :       2003 :                         mean[d] /= (double) N;
     420                 :            :                 }
     421                 :            :                 
     422                 :            :                 // Subtract data mean
     423         [ +  + ]:      51051 :                 for(int n = 0; n < N; n++) {
     424         [ +  + ]:     150200 :                         for(int d = 0; d < D; d++) {
     425                 :     100150 :                                 X[n * D + d] -= mean[d];
     426                 :            :                         }
     427                 :            :                 }
     428                 :       1001 :                 free(mean); mean = NULL;
     429                 :       1001 :         }
     430                 :            : 
     431                 :          0 :         void computeGaussianPerplexity(double* X, int N, int D, double* P, double perplexity)
     432                 :            :         {
     433                 :            :                 // Compute the squared Euclidean distance matrix
     434                 :          0 :                 double* DD = (double*) malloc(N * N * sizeof(double));
     435         [ #  # ]:          0 :                 if(DD == NULL) { printf("Memory allocation failed!\n"); exit(1); }
     436                 :          0 :                 computeSquaredEuclideanDistance(X, N, D, DD);
     437                 :            :                 
     438                 :            :                 // Compute the Gaussian kernel row by row
     439         [ #  # ]:          0 :                 for(int n = 0; n < N; n++) {
     440                 :            :                         
     441                 :            :                         // Initialize some variables
     442                 :          0 :                         bool found = false;
     443                 :          0 :                         double beta = 1.0;
     444                 :          0 :                         double min_beta = -DBL_MAX;
     445                 :          0 :                         double max_beta =  DBL_MAX;
     446                 :          0 :                         double tol = 1e-5;
     447                 :            :                         double sum_P;
     448                 :            :                         
     449                 :            :                         // Iterate until we found a good perplexity
     450                 :          0 :                         int iter = 0;
     451 [ #  # ][ #  # ]:          0 :                         while(!found && iter < 200) {
                 [ #  # ]
     452                 :            :                                 
     453                 :            :                                 // Compute Gaussian kernel row
     454         [ #  # ]:          0 :                                 for(int m = 0; m < N; m++) P[n * N + m] = exp(-beta * DD[n * N + m]);
     455                 :          0 :                                 P[n * N + n] = DBL_MIN;
     456                 :            :                                 
     457                 :            :                                 // Compute entropy of current row
     458                 :          0 :                                 sum_P = DBL_MIN;
     459         [ #  # ]:          0 :                                 for(int m = 0; m < N; m++) sum_P += P[n * N + m];
     460                 :          0 :                                 double H = 0.0;
     461         [ #  # ]:          0 :                                 for(int m = 0; m < N; m++) H += beta * (DD[n * N + m] * P[n * N + m]);
     462                 :          0 :                                 H = (H / sum_P) + log(sum_P);
     463                 :            :                                 
     464                 :            :                                 // Evaluate whether the entropy is within the tolerance level
     465                 :          0 :                                 double Hdiff = H - log(perplexity);
     466 [ #  # ][ #  # ]:          0 :                                 if(Hdiff < tol && -Hdiff < tol) {
     467                 :          0 :                                         found = true;
     468                 :            :                                 }
     469                 :            :                                 else {
     470         [ #  # ]:          0 :                                         if(Hdiff > 0) {
     471                 :          0 :                                                 min_beta = beta;
     472 [ #  # ][ #  # ]:          0 :                                                 if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
     473                 :          0 :                                                         beta *= 2.0;
     474                 :            :                                                 else
     475                 :          0 :                                                         beta = (beta + max_beta) / 2.0;
     476                 :            :                                         }
     477                 :            :                                         else {
     478                 :          0 :                                                 max_beta = beta;
     479 [ #  # ][ #  # ]:          0 :                                                 if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
     480                 :          0 :                                                         beta /= 2.0;
     481                 :            :                                                 else
     482                 :          0 :                                                         beta = (beta + min_beta) / 2.0;
     483                 :            :                                         }
     484                 :            :                                 }
     485                 :            :                                 
     486                 :            :                                 // Update iteration counter
     487                 :          0 :                                 iter++;
     488                 :            :                         }
     489                 :            :                         
     490                 :            :                         // Row normalize P
     491         [ #  # ]:          0 :                         for(int m = 0; m < N; m++) P[n * N + m] /= sum_P;
     492                 :            :                 }
     493                 :            :                 
     494                 :            :                 // Clean up memory
     495                 :          0 :                 free(DD); DD = NULL;
     496                 :          0 :         }
     497                 :            : 
     498                 :          1 :         void computeGaussianPerplexity(double* X, int N, int D, int** _row_P, int** _col_P, double** _val_P, double perplexity, int K)
     499                 :            :         { 
     500         [ -  + ]:          1 :                 if(perplexity > K) printf("Perplexity should be lower than K!\n");
     501                 :            :                 
     502                 :            :                 // Allocate the memory we need
     503                 :          1 :                 *_row_P = (int*)    malloc((N + 1) * sizeof(int));
     504                 :          1 :                 *_col_P = (int*)    calloc(N * K, sizeof(int));
     505                 :          1 :                 *_val_P = (double*) calloc(N * K, sizeof(double));
     506 [ +  - ][ +  - ]:          1 :                 if(*_row_P == NULL || *_col_P == NULL || *_val_P == NULL) { printf("Memory allocation failed!\n"); exit(1); }
                 [ -  + ]
     507                 :          1 :                 int* row_P = *_row_P;
     508                 :          1 :                 int* col_P = *_col_P;
     509                 :          1 :                 double* val_P = *_val_P;
     510                 :          1 :                 double* cur_P = (double*) malloc((N - 1) * sizeof(double));
     511         [ -  + ]:          1 :                 if(cur_P == NULL) { printf("Memory allocation failed!\n"); exit(1); }
     512                 :          1 :                 row_P[0] = 0;
     513         [ +  + ]:         51 :                 for(int n = 0; n < N; n++) row_P[n + 1] = row_P[n] + K;    
     514                 :            :                 
     515                 :            :                 // Build ball tree on data set
     516         [ +  - ]:          1 :                 VpTree<DataPoint, euclidean_distance>* tree = new VpTree<DataPoint, euclidean_distance>();
     517         [ +  - ]:          1 :                 std::vector<DataPoint> obj_X(N, DataPoint(D, -1, X));
     518 [ +  - ][ +  + ]:         51 :                 for(int n = 0; n < N; n++) obj_X[n] = DataPoint(D, n, X + n * D);
     519         [ +  - ]:          1 :                 tree->create(obj_X);
     520                 :            :                 
     521                 :            :                 // Loop over all points to find nearest neighbors
     522                 :            :                 //printf("Building tree...\n");
     523         [ +  - ]:          1 :                 std::vector<DataPoint> indices;
     524         [ +  - ]:          1 :                 std::vector<double> distances;
     525         [ +  + ]:         51 :                 for(int n = 0; n < N; n++) {
     526                 :            :                         
     527                 :            :                         //if(n % 10000 == 0) printf(" - point %d of %d\n", n, N);
     528                 :            :                         
     529                 :            :                         // Find nearest neighbors
     530         [ +  - ]:         50 :                         indices.clear();
     531         [ +  - ]:         50 :                         distances.clear();
     532         [ +  - ]:         50 :                         tree->search(obj_X[n], K + 1, &indices, &distances);
     533                 :            :                         
     534                 :            :                         // Initialize some variables for binary search
     535                 :         50 :                         bool found = false;
     536                 :         50 :                         double beta = 1.0;
     537                 :         50 :                         double min_beta = -DBL_MAX;
     538                 :         50 :                         double max_beta =  DBL_MAX;
     539                 :         50 :                         double tol = 1e-5;
     540                 :            :                         
     541                 :            :                         // Iterate until we found a good perplexity
     542                 :         50 :                         int iter = 0; double sum_P;
     543 [ +  + ][ +  - ]:       1037 :                         while(!found && iter < 200) {
                 [ +  + ]
     544                 :            :                                 
     545                 :            :                                 // Compute Gaussian kernel row
     546         [ +  + ]:      30597 :                                 for(int m = 0; m < K; m++) cur_P[m] = exp(-beta * distances[m + 1]);
     547                 :            :                                 
     548                 :            :                                 // Compute entropy of current row
     549                 :        987 :                                 sum_P = DBL_MIN;
     550         [ +  + ]:      30597 :                                 for(int m = 0; m < K; m++) sum_P += cur_P[m];
     551                 :        987 :                                 double H = .0;
     552         [ +  + ]:      30597 :                                 for(int m = 0; m < K; m++) H += beta * (distances[m + 1] * cur_P[m]);
     553                 :        987 :                                 H = (H / sum_P) + log(sum_P);
     554                 :            :                                 
     555                 :            :                                 // Evaluate whether the entropy is within the tolerance level
     556                 :        987 :                                 double Hdiff = H - log(perplexity);
     557 [ +  + ][ +  + ]:        987 :                                 if(Hdiff < tol && -Hdiff < tol) {
     558                 :         50 :                                         found = true;
     559                 :            :                                 }
     560                 :            :                                 else {
     561         [ +  + ]:        937 :                                         if(Hdiff > 0) {
     562                 :        565 :                                                 min_beta = beta;
     563 [ +  + ][ -  + ]:        565 :                                                 if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
     564                 :        255 :                                                         beta *= 2.0;
     565                 :            :                                                 else
     566                 :        565 :                                                         beta = (beta + max_beta) / 2.0;
     567                 :            :                                         }
     568                 :            :                                         else {
     569                 :        372 :                                                 max_beta = beta;
     570 [ +  - ][ -  + ]:        372 :                                                 if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
     571                 :          0 :                                                         beta /= 2.0;
     572                 :            :                                                 else
     573                 :        372 :                                                         beta = (beta + min_beta) / 2.0;
     574                 :            :                                         }
     575                 :            :                                 }
     576                 :            :                                 
     577                 :            :                                 // Update iteration counter
     578                 :        987 :                                 iter++;
     579                 :            :                         }
     580                 :            :                         
     581                 :            :                         // Row-normalize current row of P and store in matrix
     582         [ +  + ]:       1550 :                         for(int m = 0; m < K; m++) cur_P[m] /= sum_P;
     583         [ +  + ]:       1550 :                         for(int m = 0; m < K; m++) {
     584                 :       1500 :                                 col_P[row_P[n] + m] = indices[m + 1].index();
     585                 :       1500 :                                 val_P[row_P[n] + m] = cur_P[m];
     586                 :            :                         }
     587                 :            :                 }
     588                 :            :                 
     589                 :            :                 // Clean up memory
     590         [ +  - ]:          1 :                 obj_X.clear();
     591                 :          1 :                 free(cur_P);
     592 [ +  - ][ +  - ]:          1 :                 delete tree;
         [ +  - ][ +  - ]
     593                 :          1 :         }
     594                 :            : 
     595                 :            :         void computeGaussianPerplexity(double* X, int N, int D, int** _row_P, int** _col_P, double** _val_P, double perplexity, double threshold)
     596                 :            :         {
     597                 :            :                 // Allocate some memory we need for computations
     598                 :            :                 double* buff  = (double*) malloc(D * sizeof(double));
     599                 :            :                 double* DD    = (double*) malloc(N * sizeof(double));
     600                 :            :                 double* cur_P = (double*) malloc(N * sizeof(double));
     601                 :            :                 if(buff == NULL || DD == NULL || cur_P == NULL) { printf("Memory allocation failed!\n"); exit(1); }
     602                 :            : 
     603                 :            :                 // Compute the Gaussian kernel row by row (to find number of elements in sparse P)
     604                 :            :                 int total_count = 0;
     605                 :            :                 for(int n = 0; n < N; n++) {
     606                 :            :                 
     607                 :            :                         // Compute the squared Euclidean distance matrix
     608                 :            :                         for(int m = 0; m < N; m++) {
     609                 :            :                                 for(int d = 0; d < D; d++) buff[d]  = X[n * D + d];
     610                 :            :                                 for(int d = 0; d < D; d++) buff[d] -= X[m * D + d];
     611                 :            :                                 DD[m] = .0;
     612                 :            :                                 for(int d = 0; d < D; d++) DD[m] += buff[d] * buff[d];
     613                 :            :                         }
     614                 :            :                    
     615                 :            :                         // Initialize some variables
     616                 :            :                         bool found = false;
     617                 :            :                         double beta = 1.0;
     618                 :            :                         double min_beta = -DBL_MAX;
     619                 :            :                         double max_beta =  DBL_MAX;
     620                 :            :                         double tol = 1e-5;
     621                 :            :                         
     622                 :            :                         // Iterate until we found a good perplexity
     623                 :            :                         int iter = 0; double sum_P;
     624                 :            :                         while(!found && iter < 200) {
     625                 :            :                                 
     626                 :            :                                 // Compute Gaussian kernel row
     627                 :            :                                 for(int m = 0; m < N; m++) cur_P[m] = exp(-beta * DD[m]);
     628                 :            :                                 cur_P[n] = DBL_MIN;
     629                 :            :                                 
     630                 :            :                                 // Compute entropy of current row
     631                 :            :                                 sum_P = DBL_MIN;
     632                 :            :                                 for(int m = 0; m < N; m++) sum_P += cur_P[m];
     633                 :            :                                 double H = 0.0;
     634                 :            :                                 for(int m = 0; m < N; m++) H += beta * (DD[m] * cur_P[m]);
     635                 :            :                                 H = (H / sum_P) + log(sum_P);
     636                 :            :                                 
     637                 :            :                                 // Evaluate whether the entropy is within the tolerance level
     638                 :            :                                 double Hdiff = H - log(perplexity);
     639                 :            :                                 if(Hdiff < tol && -Hdiff < tol) {
     640                 :            :                                         found = true;
     641                 :            :                                 }
     642                 :            :                                 else {
     643                 :            :                                         if(Hdiff > 0) {
     644                 :            :                                                 min_beta = beta;
     645                 :            :                                                 if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
     646                 :            :                                                         beta *= 2.0;
     647                 :            :                                                 else
     648                 :            :                                                         beta = (beta + max_beta) / 2.0;
     649                 :            :                                         }
     650                 :            :                                         else {
     651                 :            :                                                 max_beta = beta;
     652                 :            :                                                 if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
     653                 :            :                                                         beta /= 2.0;
     654                 :            :                                                 else
     655                 :            :                                                         beta = (beta + min_beta) / 2.0;
     656                 :            :                                         }
     657                 :            :                                 }
     658                 :            :                                 
     659                 :            :                                 // Update iteration counter
     660                 :            :                                 iter++;
     661                 :            :                         }
     662                 :            :                         
     663                 :            :                         // Row-normalize and threshold current row of P
     664                 :            :                         for(int m = 0; m < N; m++) cur_P[m] /= sum_P;
     665                 :            :                         for(int m = 0; m < N; m++) {
     666                 :            :                                 if(cur_P[m] > threshold / (double) N) total_count++;
     667                 :            :                         }
     668                 :            :                 }
     669                 :            :                 
     670                 :            :                 // Allocate the memory we need
     671                 :            :                 *_row_P = (int*)    malloc((N + 1)     * sizeof(int));
     672                 :            :                 *_col_P = (int*)    malloc(total_count * sizeof(int));
     673                 :            :                 *_val_P = (double*) malloc(total_count * sizeof(double));
     674                 :            :                 int* row_P = *_row_P;
     675                 :            :                 int* col_P = *_col_P;
     676                 :            :                 double* val_P = *_val_P;
     677                 :            :                 if(row_P == NULL || col_P == NULL || val_P == NULL) { printf("Memory allocation failed!\n"); exit(1); }
     678                 :            :                 row_P[0] = 0;
     679                 :            :                 
     680                 :            :                 // Compute the Gaussian kernel row by row (this time, store the results)
     681                 :            :                 int count = 0;
     682                 :            :                 for(int n = 0; n < N; n++) {
     683                 :            :                         
     684                 :            :                         // Compute the squared Euclidean distance matrix
     685                 :            :                         for(int m = 0; m < N; m++) {
     686                 :            :                                 for(int d = 0; d < D; d++) buff[d]  = X[n * D + d];
     687                 :            :                                 for(int d = 0; d < D; d++) buff[d] -= X[m * D + d];
     688                 :            :                                 DD[m] = .0;
     689                 :            :                                 for(int d = 0; d < D; d++) DD[m] += buff[d] * buff[d];
     690                 :            :                         }
     691                 :            :                         
     692                 :            :                         // Initialize some variables
     693                 :            :                         bool found = false;
     694                 :            :                         double beta = 1.0;
     695                 :            :                         double min_beta = -DBL_MAX;
     696                 :            :                         double max_beta =  DBL_MAX;
     697                 :            :                         double tol = 1e-5;
     698                 :            :                         
     699                 :            :                         // Iterate until we found a good perplexity
     700                 :            :                         int iter = 0; double sum_P;
     701                 :            :                         while(!found && iter < 200) {
     702                 :            :                                 
     703                 :            :                                 // Compute Gaussian kernel row
     704                 :            :                                 for(int m = 0; m < N; m++) cur_P[m] = exp(-beta * DD[m]);
     705                 :            :                                 cur_P[n] = DBL_MIN;
     706                 :            :                                 
     707                 :            :                                 // Compute entropy of current row
     708                 :            :                                 sum_P = DBL_MIN;
     709                 :            :                                 for(int m = 0; m < N; m++) sum_P += cur_P[m];
     710                 :            :                                 double H = 0.0;
     711                 :            :                                 for(int m = 0; m < N; m++) H += beta * (DD[m] * cur_P[m]);
     712                 :            :                                 H = (H / sum_P) + log(sum_P);
     713                 :            :                                 
     714                 :            :                                 // Evaluate whether the entropy is within the tolerance level
     715                 :            :                                 double Hdiff = H - log(perplexity);
     716                 :            :                                 if(Hdiff < tol && -Hdiff < tol) {
     717                 :            :                                         found = true;
     718                 :            :                                 }
     719                 :            :                                 else {
     720                 :            :                                         if(Hdiff > 0) {
     721                 :            :                                                 min_beta = beta;
     722                 :            :                                                 if(max_beta == DBL_MAX || max_beta == -DBL_MAX)
     723                 :            :                                                         beta *= 2.0;
     724                 :            :                                                 else
     725                 :            :                                                         beta = (beta + max_beta) / 2.0;
     726                 :            :                                         }
     727                 :            :                                         else {
     728                 :            :                                                 max_beta = beta;
     729                 :            :                                                 if(min_beta == -DBL_MAX || min_beta == DBL_MAX)
     730                 :            :                                                         beta /= 2.0;
     731                 :            :                                                 else
     732                 :            :                                                         beta = (beta + min_beta) / 2.0;
     733                 :            :                                         }
     734                 :            :                                 }
     735                 :            :                                 
     736                 :            :                                 // Update iteration counter
     737                 :            :                                 iter++;
     738                 :            :                         }
     739                 :            :                         
     740                 :            :                         // Row-normalize and threshold current row of P
     741                 :            :                         for(int m = 0; m < N; m++) cur_P[m] /= sum_P;
     742                 :            :                         for(int m = 0; m < N; m++) {
     743                 :            :                                 if(cur_P[m] > threshold / (double) N) {
     744                 :            :                                         col_P[count] = m;
     745                 :            :                                         val_P[count] = cur_P[m];
     746                 :            :                                         count++;
     747                 :            :                                 }
     748                 :            :                         }
     749                 :            :                         row_P[n + 1] = count;
     750                 :            :                 }
     751                 :            :                 
     752                 :            :                 // Clean up memory
     753                 :            :                 free(DD);    DD = NULL;
     754                 :            :                 free(buff);  buff = NULL;
     755                 :            :                 free(cur_P); cur_P = NULL;
     756                 :            :         }
     757                 :            : 
     758                 :          0 :         void computeSquaredEuclideanDistance(double* X, int N, int D, double* DD)
     759                 :            :         {
     760                 :          0 :                 double* dataSums = (double*) calloc(N, sizeof(double));
     761         [ #  # ]:          0 :                 if(dataSums == NULL) { printf("Memory allocation failed!\n"); exit(1); }
     762         [ #  # ]:          0 :                 for(int n = 0; n < N; n++) {
     763         [ #  # ]:          0 :                         for(int d = 0; d < D; d++) {
     764                 :          0 :                                 dataSums[n] += (X[n * D + d] * X[n * D + d]);
     765                 :            :                         }
     766                 :            :                 }
     767         [ #  # ]:          0 :                 for(int n = 0; n < N; n++) {
     768         [ #  # ]:          0 :                         for(int m = 0; m < N; m++) {
     769                 :          0 :                                 DD[n * N + m] = dataSums[n] + dataSums[m];
     770                 :            :                         }
     771                 :            :                 }
     772                 :          0 :                 Eigen::Map<Eigen::MatrixXd> DD_map(DD,N,N);
     773                 :          0 :                 Eigen::Map<Eigen::MatrixXd> X_map(X,D,N);
     774 [ #  # ][ #  # ]:          0 :                 DD_map.noalias() = -2.0*X_map.transpose()*X_map;
     775                 :            : 
     776                 :            :                 //cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans, N, N, D, -2.0, X, D, X, D, 1.0, DD, N);
     777                 :          0 :                 free(dataSums); dataSums = NULL;
     778                 :          0 :         }
     779                 :            : 
     780                 :            : };
     781                 :            : 
     782                 :            : }
     783                 :            : 
     784                 :            : #endif

Generated by: LCOV version 1.9