diff --git a/Makefile b/Makefile index be317d5..b86f3a4 100644 --- a/Makefile +++ b/Makefile @@ -1,8 +1,10 @@ -build: +build: tsne/bh_sne.pyx \ + tsne/bh_sne_src/tsne.cpp \ + tsne/bh_sne_src/sptree.cpp \ + $(wildcard tsne/bh_sne_src/*.h) python setup.py build_ext --inplace -install: - python setup.py build_ext --inplace +install: build python setup.py install sdist: diff --git a/README.md b/README.md index c0d577e..169dad0 100644 --- a/README.md +++ b/README.md @@ -23,7 +23,6 @@ Requirements * [numpy](numpy.scipy.org) > =1.7.1 * [scipy](http://www.scipy.org/) >= 0.12.0 * [cython](cython.org) >= 0.19.1 -* [cblas](http://www.netlib.org/blas/) or [openblas](https://github.com/xianyi/OpenBLAS). Tested version is v0.2.5 and v0.2.6 (not necessary for OSX). [Anaconda](http://continuum.io/downloads) is recommended. diff --git a/setup.py b/setup.py index 7822cf6..50ad878 100644 --- a/setup.py +++ b/setup.py @@ -35,8 +35,8 @@ ext_modules = [Extension(name='bh_sne', sources=['tsne/bh_sne_src/sptree.cpp', 'tsne/bh_sne_src/tsne.cpp', 'tsne/bh_sne.pyx'], include_dirs=[numpy.get_include(), 'tsne/bh_sne_src/'], - extra_compile_args=extra_compile_args + ['-ffast-math', '-O3'], - extra_link_args=['-Wl,-framework', '-Wl,Accelerate', '-lcblas'], + extra_compile_args=extra_compile_args + ['-ffast-math', '-O3', '-std=c++11'], + extra_link_args=['-Wl,-framework', '-Wl,Accelerate'], language='c++')] else: @@ -44,18 +44,10 @@ ext_modules = [Extension(name='bh_sne', sources=['tsne/bh_sne_src/sptree.cpp', 'tsne/bh_sne_src/tsne.cpp', 'tsne/bh_sne.pyx'], - include_dirs=[numpy.get_include(), '/usr/local/include', 'tsne/bh_sne_src/'], - library_dirs=['/usr/local/lib', '/usr/lib64/atlas'], - extra_compile_args=['-msse2', '-O3', '-fPIC', '-w', '-ffast-math'], - extra_link_args=['-Wl,-Bstatic', '-lcblas', '-Wl,-Bdynamic'], - language='c++'), - - Extension(name='bh_sne_3d', - sources=['tsne/bh_sne_src/sptree.cpp', 'tsne/bh_sne_src/tsne.cpp', 'tsne/bh_sne_3d.pyx'], - include_dirs=[numpy.get_include(), '/usr/local/include', 'tsne/bh_sne_src/'], - library_dirs=['/usr/local/lib', '/usr/lib64/atlas'], - extra_compile_args=['-msse2', '-O3', '-fPIC', '-w', '-ffast-math', '-DTSNE3D'], - extra_link_args=['-Wl,-Bstatic', '-lcblas', '-Wl,-Bdynamic'], + include_dirs=[numpy.get_include(), 'tsne/bh_sne_src/'], + extra_compile_args=['-msse3', '-O3', '-fPIC', '-w', '-ffast-math', '-std=c++11', + '-ffunction-sections', '-flto', '-mtune=native'], + extra_link_args=['-O3', '-Wl,--gc-sections', '-flto', '-mtune=native'], language='c++')] ext_modules = cythonize(ext_modules) diff --git a/tsne/__init__.py b/tsne/__init__.py index c5b9d3c..4cc3b0f 100644 --- a/tsne/__init__.py +++ b/tsne/__init__.py @@ -3,8 +3,7 @@ import numpy as np import scipy.linalg as la import sys -from bh_sne import BH_SNE -from bh_sne_3d import BH_SNE_3D +from bh_sne import BH_SNE as tsne def bh_sne(data, pca_d=None, d=2, perplexity=30., theta=0.5, random_state=None, copy_data=False, init=None, @@ -79,13 +78,6 @@ def bh_sne(data, pca_d=None, d=2, perplexity=30., theta=0.5, if mom_switch_iter is None: mom_switch_iter = 250 - if d == 2: - tsne = BH_SNE() - elif d == 3: - tsne = BH_SNE_3D() - else: - raise Exception("TSNE dimensions must be 2 or 3") - Y = tsne.run(X, N, X.shape[1], d, perplexity, theta, seed, init=init, use_init=use_init, max_iter=max_iter, stop_lying_iter=stop_lying_iter, mom_switch_iter=mom_switch_iter) return Y diff --git a/tsne/bh_sne.pyx b/tsne/bh_sne.pyx index e45dd13..77e4832 100644 --- a/tsne/bh_sne.pyx +++ b/tsne/bh_sne.pyx @@ -4,25 +4,35 @@ cimport numpy as np cimport cython from libcpp cimport bool -cdef extern from "tsne.h": - cdef cppclass TSNE: - TSNE() - void run(double* X, int N, int D, double* Y, int no_dims, double perplexity, double theta, int rand_seed, bool skip_random_init, double *init, bool use_init, int max_iter, int stop_lying_iter, int mom_switch_iter) +cdef extern from "tsne.h" namespace "TSNE": + void c_run "TSNE::run" (double* X, int N, int D, double* Y, int no_dims, + double perplexity, double theta, + int rand_seed, bool skip_random_init, + double *init, bool use_init, + int max_iter, int stop_lying_iter, int mom_switch_iter) nogil cdef class BH_SNE: - cdef TSNE* thisptr # hold a C++ instance - - def __cinit__(self): - self.thisptr = new TSNE() - - def __dealloc__(self): - del self.thisptr - @cython.boundscheck(False) @cython.wraparound(False) - def run(self, X, N, D, d, perplexity, theta, seed, init, use_init, max_iter, stop_lying_iter, mom_switch_iter): - cdef np.ndarray[np.float64_t, ndim=2, mode='c'] _X = np.ascontiguousarray(X) - cdef np.ndarray[np.float64_t, ndim=2, mode='c'] _init = np.ascontiguousarray(init) - cdef np.ndarray[np.float64_t, ndim=2, mode='c'] Y = np.zeros((N, d), dtype=np.float64) - self.thisptr.run(&_X[0,0], N, D, &Y[0,0], d, perplexity, theta, seed, False, &_init[0,0], use_init, max_iter, stop_lying_iter, mom_switch_iter) + @staticmethod + def run(X, int N, int D, int d, + double perplexity, double theta, + int seed, init, bool use_init, + int max_iter, int stop_lying_iter, int mom_switch_iter): + cdef np.ndarray[np.float64_t, ndim=2, mode='c'] _X = np.ascontiguousarray( + X, + dtype=np.float64) + cdef np.ndarray[np.float64_t, ndim=2, mode='c'] _init = np.ascontiguousarray( + init, + dtype=np.float64) + cdef np.ndarray[np.float64_t, ndim=2, mode='c'] Y = np.zeros( + (N, d), + dtype=np.float64, + order='C') + assert(N, X.shape[1]) + with nogil: + c_run(&_X[0,0], N, D, &Y[0,0], d, + perplexity, theta, + seed, False, &_init[0,0], use_init, + max_iter, stop_lying_iter, mom_switch_iter) return Y diff --git a/tsne/bh_sne_3d.pyx b/tsne/bh_sne_3d.pyx deleted file mode 100644 index 315e107..0000000 --- a/tsne/bh_sne_3d.pyx +++ /dev/null @@ -1,28 +0,0 @@ -# distutils: language = c++ -import numpy as np -cimport numpy as np -cimport cython -from libcpp cimport bool - -cdef extern from "tsne.h": - cdef cppclass TSNE: - TSNE() - void run(double* X, int N, int D, double* Y, int no_dims, double perplexity, double theta, int rand_seed, bool skip_random_init, double *init, bool use_init, int max_iter, int stop_lying_iter, int mom_switch_iter) - -cdef class BH_SNE_3D: - cdef TSNE* thisptr # hold a C++ instance - - def __cinit__(self): - self.thisptr = new TSNE() - - def __dealloc__(self): - del self.thisptr - - @cython.boundscheck(False) - @cython.wraparound(False) - def run(self, X, N, D, d, perplexity, theta, seed, init, use_init, max_iter, stop_lying_iter, mom_switch_iter): - cdef np.ndarray[np.float64_t, ndim=2, mode='c'] _X = np.ascontiguousarray(X) - cdef np.ndarray[np.float64_t, ndim=2, mode='c'] _init = np.ascontiguousarray(init) - cdef np.ndarray[np.float64_t, ndim=2, mode='c'] Y = np.zeros((N, d), dtype=np.float64) - self.thisptr.run(&_X[0,0], N, D, &Y[0,0], d, perplexity, theta, seed, False, &_init[0,0], use_init, max_iter, stop_lying_iter, mom_switch_iter) - return Y diff --git a/tsne/bh_sne_src/.clang-format b/tsne/bh_sne_src/.clang-format new file mode 100644 index 0000000..9ebca72 --- /dev/null +++ b/tsne/bh_sne_src/.clang-format @@ -0,0 +1,3 @@ +BasedOnStyle: Google +IndentWidth: 4 +AllowShortIfStatementsOnASingleLine: false diff --git a/tsne/bh_sne_src/Makefile b/tsne/bh_sne_src/Makefile index 8c8a1c8..42b4bcf 100644 --- a/tsne/bh_sne_src/Makefile +++ b/tsne/bh_sne_src/Makefile @@ -2,25 +2,28 @@ #CFLAGS = -march=haswell -ffast-math -O3 -Rpass=loop-vectorize -Rpass-missed=loop-vectorize -Rpass-analysis=loop-vectorize #CFLAGS = -march=haswell -ffast-math -O3 -CXX = g++ -CFLAGS = -ffast-math -O3 +CXX ?= g++ +CFLAGS += -std=c++11 -ffast-math -O3 -flto -all: bh_tsne bh_tsne_3d +CFLAGS += -ffunction-sections +LDFLAGS += -Wl,--gc-sections -bh_tsne: tsne.o sptree.o - $(CXX) $(CFLAGS) tsne.o sptree.o -o bh_tsne +all: bh_tsne -bh_tsne_3d: tsne_3d.o sptree.o - $(CXX) $(CFLAGS) tsne_3d.o sptree.o -o bh_tsne_3d +bh_tsne: main.o io.o tsne.o sptree.o + $(CXX) $(CFLAGS) $(LDFLAGS) $^ -o $@ sptree.o: sptree.cpp sptree.h - $(CXX) $(CFLAGS) -c sptree.cpp + $(CXX) $(CFLAGS) -c $< -o $@ tsne.o: tsne.cpp tsne.h sptree.h vptree.h - $(CXX) $(CFLAGS) -c tsne.cpp + $(CXX) $(CFLAGS) -c $< -o $@ -tsne_3d.o: tsne.cpp tsne.h sptree.h vptree.h - $(CXX) $(CFLAGS) -DTSNE3D -c tsne.cpp +main.o: main.cpp tsne.h sptree.h vptree.h + $(CXX) $(CFLAGS) -c $< -o $@ + +io.o: io.cpp tsne.h + $(CXX) $(CFLAGS) -c $< -o $@ clean: - rm -Rf *.o bh_tsne bh_tsne_3d + rm -Rf *.o bh_tsne diff --git a/tsne/bh_sne_src/io.cpp b/tsne/bh_sne_src/io.cpp new file mode 100644 index 0000000..3944c19 --- /dev/null +++ b/tsne/bh_sne_src/io.cpp @@ -0,0 +1,103 @@ +/* + * + * Copyright (c) 2014, Laurens van der Maaten (Delft University of Technology) + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * 3. All advertising materials mentioning features or use of this software + * must display the following acknowledgement: + * This product includes software developed by the Delft University of Technology. + * 4. Neither the name of the Delft University of Technology nor the names of + * its contributors may be used to endorse or promote products derived from + * this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY LAURENS VAN DER MAATEN ''AS IS'' AND ANY EXPRESS + * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO + * EVENT SHALL LAURENS VAN DER MAATEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR + * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * 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. + * + */ + +#include +#include +#include + +#include "tsne.h" + +namespace TSNE { + +// Function that loads data from a t-SNE file +// Note: this function does a malloc that should be freed elsewhere +bool load_data(const char* dat_file, double** data, int* n, int* d, int* no_dims, double* theta, double* perplexity, int* rand_seed, int* max_iter) { + + // Open file, read first 2 integers, allocate memory, and read the data + FILE *h; + if((h = fopen(dat_file, "r+b")) == NULL) { + fprintf(stderr,"Error: could not open data file.\n"); + return false; + } + fread(n, sizeof(int), 1, h); // number of datapoints + fread(d, sizeof(int), 1, h); // original dimensionality + fread(theta, sizeof(double), 1, h); // gradient accuracy + fread(perplexity, sizeof(double), 1, h); // perplexity + fread(no_dims, sizeof(int), 1, h); // output dimensionality + fread(max_iter, sizeof(int),1,h); // maximum number of iterations + *data = (double*) malloc(*d * *n * sizeof(double)); + if(*data == NULL) { fprintf(stderr,"Memory allocation failed!\n"); exit(1); } + fread(*data, sizeof(double), *n * *d, h); // the data + if(!feof(h)) fread(rand_seed, sizeof(int), 1, h); // random seed + fclose(h); + fprintf(stderr,"Read the %i x %i data matrix successfully!\n", *n, *d); + return true; +} + +// Function that saves map to a t-SNE file +void save_data(const char* res_file, double* data, int* landmarks, double* costs, int n, int d) { + + // Open file, write first 2 integers and then the data + FILE *h; + if((h = fopen(res_file, "w+b")) == NULL) { + fprintf(stderr,"Error: could not open data file.\n"); + return; + } + fwrite(&n, sizeof(int), 1, h); + fwrite(&d, sizeof(int), 1, h); + fwrite(data, sizeof(double), n * d, h); + fwrite(landmarks, sizeof(int), n, h); + fwrite(costs, sizeof(double), n, h); + fclose(h); + fprintf(stderr,"Wrote the %i x %i data matrix successfully!\n", n, d); +} + +void save_csv(const char* csv_file, double* Y, int N, int D) { + std::ofstream csv(csv_file); + + for (int d = 0; d < D; d++) { + csv << "TSNE" << d+1 << ","; + } + csv << "\n"; + + for (int n = 0; n < N; n++) { + int row_offset = n * D; + for (int d = 0; d < D; d++) { + csv << Y[row_offset + d] << ","; + } + csv << "\n"; + } + + csv.close(); +} + +} // namespace TSNE diff --git a/tsne/bh_sne_src/main.cpp b/tsne/bh_sne_src/main.cpp new file mode 100644 index 0000000..86d8205 --- /dev/null +++ b/tsne/bh_sne_src/main.cpp @@ -0,0 +1,80 @@ +/* + * + * Copyright (c) 2014, Laurens van der Maaten (Delft University of Technology) + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * 3. All advertising materials mentioning features or use of this software + * must display the following acknowledgement: + * This product includes software developed by the Delft University of Technology. + * 4. Neither the name of the Delft University of Technology nor the names of + * its contributors may be used to endorse or promote products derived from + * this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY LAURENS VAN DER MAATEN ''AS IS'' AND ANY EXPRESS + * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO + * EVENT SHALL LAURENS VAN DER MAATEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR + * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * 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. + * + */ + + +#include "tsne.h" + +#include +#include +#include + +using std::fprintf; +using std::free; +using std::vector; + +// Function that runs the Barnes-Hut implementation of t-SNE +int main(int argc, char *argv[]) { + + // load input and output + const char *dat_file = "data.dat"; + const char *res_file = "result.dat"; + if (argc > 1) { + dat_file = argv[1]; + res_file = argv[2]; + } + + // Define some variables + int origN, N, D, no_dims, max_iter, *landmarks; + double perc_landmarks; + double perplexity, theta, *data; + int rand_seed = -1; + + // Read the parameters and the dataset + if(TSNE::load_data(dat_file, &data, &origN, &D, &no_dims, &theta, &perplexity, &rand_seed, &max_iter)) { + + // Make dummy landmarks + N = origN; + + // Now fire up the SNE implementation + vector Y(N * no_dims); + TSNE::run(data, N, D, Y.data(), no_dims, perplexity, theta, rand_seed, false, NULL, false, max_iter); + + // Save the results + vector landmarks(N); + for(int n = 0; n < N; n++) landmarks[n] = n; + vector costs(N); + TSNE::save_data(res_file, Y.data(), landmarks.data(), costs.data(), N, no_dims); + + // Clean up the memory + free(data); data = NULL; + } +} diff --git a/tsne/bh_sne_src/sptree.cpp b/tsne/bh_sne_src/sptree.cpp index 2d187c6..40f6f1a 100644 --- a/tsne/bh_sne_src/sptree.cpp +++ b/tsne/bh_sne_src/sptree.cpp @@ -30,421 +30,426 @@ * */ -#include -#include -#include -#include -#include #include "sptree.h" +#include // for max, max_element +#include // for DBL_MAX +#include +#include // for fprintf, stderr, size_t +#include // for unique_ptr +using std::abs; +using std::array; +using std::max; +using std::max_element; +using std::unique_ptr; +using std::vector; -// Constructs cell -template -Cell::Cell() { -} +namespace TSNE { +namespace _sptree_internal { -template -Cell::Cell(double* inp_corner, double* inp_width) { - for(int d = 0; d < NDims; d++) setCorner(d, inp_corner[d]); - for(int d = 0; d < NDims; d++) setWidth( d, inp_width[d]); -} +template +Cell::Cell(const typename Cell::point_t& p) : corner(p) {} -// Destructs cell -template -Cell::~Cell() { -} - -template -double Cell::getCorner(unsigned int d) { +template +double Cell::getCorner(unsigned int d) const { return corner[d]; } -template -double Cell::getWidth(unsigned int d) { - return width[d]; +template +void Cell::setCorner(const typename Cell::point_t& val) { + corner = val; } -template +template void Cell::setCorner(unsigned int d, double val) { corner[d] = val; } -template -void Cell::setWidth(unsigned int d, double val) { - width[d] = val; -} - // Checks whether a point lies in a cell -template -bool Cell::containsPoint(double point[]) -{ - for(int d = 0; d < NDims; d++) { - if(corner[d] - width[d] > point[d]) return false; - if(corner[d] + width[d] < point[d]) return false; +template +bool Cell::containsPoint( + const double* point, const typename Cell::point_t& width) const { + for (int d = 0; d < NDims; ++d) { + if (abs(corner[d] - point[d]) > width[d]) + return false; } return true; } - -// Default constructor for SPTree -- build tree, too! -template -SPTree::SPTree(double* inp_data, unsigned int N) -{ - - // Compute mean, width, and height of current map (boundaries of SPTree) - int nD = 0; - double* mean_Y = (double*) calloc(NDims, sizeof(double)); - double* min_Y = (double*) malloc(NDims * sizeof(double)); - double* max_Y = (double*) malloc(NDims * sizeof(double)); - - for(unsigned int d = 0; d < NDims; d++) { - min_Y[d] = DBL_MAX; - max_Y[d] = -DBL_MAX; - } - - for(unsigned int n = 0; n < N; n++) { - for(unsigned int d = 0; d < NDims; d++) { - mean_Y[d] += inp_data[n * NDims + d]; - if(inp_data[nD + d] < min_Y[d]) min_Y[d] = inp_data[nD + d]; - if(inp_data[nD + d] > max_Y[d]) max_Y[d] = inp_data[nD + d]; - } - nD += NDims; - } - - for(int d = 0; d < NDims; d++) mean_Y[d] /= (double) N; - - // Construct SPTree - double* width = (double*) malloc(NDims * sizeof(double)); - for(int d = 0; d < NDims; d++) width[d] = fmax(max_Y[d] - mean_Y[d], mean_Y[d] - min_Y[d]) + 1e-5; - init(NULL, inp_data, mean_Y, width); - fill(N); - - // Clean up memory - free(mean_Y); - free(max_Y); - free(min_Y); - free(width); +// Constructor for SPTreeNode. +template +SPTreeNode::SPTreeNode( + const typename SPTreeNode::point_t& inp_corner) + : boundary(inp_corner), cum_size(0) { + center_of_mass.fill(0.0); } - - -// Constructor for SPTree with particular size and parent -- build the tree, too! -template -SPTree::SPTree(double* inp_data, unsigned int N, double* inp_corner, double* inp_width) -{ - init(NULL, inp_data, inp_corner, inp_width); - fill(N); +// Constructor for SPTreeNode. +template +SPTreeNode::SPTreeNode() : cum_size(0) { + center_of_mass.fill(0.0); } - -// Constructor for SPTree with particular size (do not fill the tree) -template -SPTree::SPTree(double* inp_data, double* inp_corner, double* inp_width) -{ - init(NULL, inp_data, inp_corner, inp_width); +// Update the corner position. +template +void SPTreeNode::setCorner( + const typename SPTreeNode::point_t& inp_corner) { + boundary.setCorner(inp_corner); } - -// Constructor for SPTree with particular size and parent (do not fill tree) -template -SPTree::SPTree(SPTree* inp_parent, double* inp_data, double* inp_corner, double* inp_width) { - init(inp_parent, inp_data, inp_corner, inp_width); +template +bool SPTreeNode::is_leaf() const { + return NDims == 0 || cum_size <= QT_NODE_CAPACITY; } - -// Constructor for SPTree with particular size and parent -- build the tree, too! -template -SPTree::SPTree(SPTree* inp_parent, double* inp_data, unsigned int N, double* inp_corner, double* inp_width) -{ - init(inp_parent, inp_data, inp_corner, inp_width); - fill(N); -} - - -// Main initialization function -template -void SPTree::init(SPTree* inp_parent, double* inp_data, double* inp_corner, double* inp_width) -{ - parent = inp_parent; - data = inp_data; - is_leaf = true; - size = 0; - cum_size = 0; - - for(unsigned int d = 0; d < NDims; d++) boundary.setCorner(d, inp_corner[d]); - for(unsigned int d = 0; d < NDims; d++) boundary.setWidth( d, inp_width[d]); - - for(unsigned int i = 0; i < no_children; i++) children[i] = NULL; - for(unsigned int d = 0; d < NDims; d++) center_of_mass[d] = .0; -} - - -// Destructor for SPTree -template -SPTree::~SPTree() -{ - for(unsigned int i = 0; i < no_children; i++) { - if(children[i] != NULL) delete children[i]; +template +unsigned int SPTreeNode::which_child(const double* point) const { + unsigned int div = 1; + unsigned int i = 0; + for (int d = 0; d < NDims; d++) { + if (boundary.getCorner(d) > point[d]) { + i |= div; + } + div *= 2; } + return i; } - -// Update the data underlying this tree -template -void SPTree::setData(double* inp_data) -{ - data = inp_data; -} - - -// Get the parent of the current tree -template -SPTree* SPTree::getParent() -{ - return parent; -} - - -// Insert a point into the SPTree -template -bool SPTree::insert(unsigned int new_index) -{ +template +bool SPTreeNode::insert(unsigned int new_index, const double* data, + vector* widths, + typename vector::size_type depth) { // Ignore objects which do not belong in this quad tree - double* point = data + new_index * NDims; - if(!boundary.containsPoint(point)) + const double* point = data + new_index * NDims; + if (!boundary.containsPoint(point, (*widths)[depth])) return false; // Online update of cumulative size and center-of-mass - cum_size++; - double mult1 = (double) (cum_size - 1) / (double) cum_size; - double mult2 = 1.0 / (double) cum_size; + unsigned int size = cum_size++; + if (size != 0) { + double mult2 = 1.0 / (double)cum_size; + double mult1 = (double)(size)*mult2; - for(unsigned int d = 0; d < NDims; d++) { - center_of_mass[d] = center_of_mass[d] * mult1 + mult2 * point[d]; + for (unsigned int d = 0; d < NDims; d++) { + center_of_mass[d] = center_of_mass[d] * mult1 + mult2 * point[d]; + } + } else { + for (unsigned int d = 0; d < NDims; d++) { + center_of_mass[d] = point[d]; + } } // If there is space in this quad tree and it is a leaf, add the object here - if(is_leaf && size < QT_NODE_CAPACITY) { + if (size < QT_NODE_CAPACITY) { index[size] = new_index; - size++; return true; - } - - // Don't add duplicates for now (this is not very nice) - bool any_duplicate = false; - for(unsigned int n = 0; n < size; n++) { - bool duplicate = true; - for(unsigned int d = 0; d < NDims; d++) { - if(point[d] != data[index[n] * NDims + d]) { duplicate = false; break; } + } else if (size == QT_NODE_CAPACITY) { + // Don't add duplicates for now (this is not very nice) + for (unsigned int n = 0; n < size; n++) { + if (__builtin_expect(index[n] == new_index, 0)) { + cum_size--; + return true; + } + bool duplicate = true; + const double* dp = data + index[n] * NDims; + for (unsigned int d = 0; d < NDims; d++) { + if (__builtin_expect(point[d] != dp[d], 1)) { + duplicate = false; + break; + } + } + if (__builtin_expect(duplicate, 0)) { + cum_size--; + return true; + } } - any_duplicate = any_duplicate | duplicate; + // We need to subdivide the current cell + subdivide(data, widths, depth); } - if(any_duplicate) return true; - - // Otherwise, we need to subdivide the current cell - if(is_leaf) subdivide(); // Find out where the point can be inserted - for(unsigned int i = 0; i < no_children; i++) { - if(children[i]->insert(new_index)) return true; - } - - // Otherwise, the point cannot be inserted (this should never happen) - return false; + auto c = which_child(point); + return (*children)[c].insert(new_index, data, widths, depth + 1); } - -// Create four children which fully divide this cell into four quads of equal area -template -void SPTree::subdivide() { +// Create four children which fully divide this cell into four quads of equal +// area +template +void SPTreeNode::subdivide(const double* data, vector* widths, + typename vector::size_type depth) { + // If nessessary, add to the width. + if (depth + 1 == widths->size()) { + // extend the list. + widths->emplace_back(); + point_t& child_widths = widths->back(); + const point_t& width = (*widths)[depth]; + for (unsigned int d = 0; d < NDims; d++) { + child_widths[d] = .5 * width[d]; + } + } + const point_t& new_width = (*widths)[depth + 1]; // Create new children - double new_corner[NDims]; - double new_width[NDims]; - for(unsigned int i = 0; i < no_children; i++) { - unsigned int div = 1; - for(unsigned int d = 0; d < NDims; d++) { - new_width[d] = .5 * boundary.getWidth(d); - if((i / div) % 2 == 1) new_corner[d] = boundary.getCorner(d) - .5 * boundary.getWidth(d); - else new_corner[d] = boundary.getCorner(d) + .5 * boundary.getWidth(d); - div *= 2; + children.reset(new array()); + auto& chi = (*children); + for (unsigned int i = 0; i < no_children; i++) { + Cell& new_corner = chi[i].boundary; + for (unsigned int d = 0; d < NDims; d++) { + if ((i >> d) % 2 == 1) + new_corner.setCorner(d, boundary.getCorner(d) - new_width[d]); + else + new_corner.setCorner(d, boundary.getCorner(d) + new_width[d]); } - children[i] = new SPTree(this, data, new_corner, new_width); } // Move existing points to correct children - for(unsigned int i = 0; i < size; i++) { - bool success = false; - for(unsigned int j = 0; j < no_children; j++) { - if(!success) success = children[j]->insert(index[i]); - } - index[i] = -1; + for (unsigned int i = 0; i < QT_NODE_CAPACITY; i++) { + auto c = which_child(data + index[i] * NDims); + chi[c].insert(index[i], data, widths, depth + 1); } - - // Empty parent node - size = 0; - is_leaf = false; } - -// Build SPTree on dataset -template -void SPTree::fill(unsigned int N) -{ - for(unsigned int i = 0; i < N; i++) insert(i); -} - - -// Checks whether the specified tree is correct -template -bool SPTree::isCorrect() -{ - for(unsigned int n = 0; n < size; n++) { - double* point = data + index[n] * NDims; - if(!boundary.containsPoint(point)) return false; - } - if(!is_leaf) { - bool correct = true; - for(int i = 0; i < no_children; i++) correct = correct && children[i]->isCorrect(); - return correct; +template +bool SPTreeNode::isCorrect( + const double* data, typename vector::const_iterator width) const { + if (is_leaf()) { + for (unsigned int n = 0; n < cum_size; n++) { + const double* point = data + index[n] * NDims; + if (!boundary.containsPoint(point, *width)) + return false; + } + } else { + ++width; + for (const auto& child : *children) { + if (!child.isCorrect(data, width)) { + return false; + } + } } - else return true; -} - - - -// Build a list of all indices in SPTree -template -void SPTree::getAllIndices(unsigned int* indices) -{ - getAllIndices(indices, 0); + return true; } - // Build a list of all indices in SPTree -template -unsigned int SPTree::getAllIndices(unsigned int* indices, unsigned int loc) -{ - +template +unsigned int SPTreeNode::getAllIndices(unsigned int* indices, + unsigned int loc) const { // Gather indices in current quadrant - for(unsigned int i = 0; i < size; i++) indices[loc + i] = index[i]; - loc += size; - - // Gather indices in children - if(!is_leaf) { - for(int i = 0; i < no_children; i++) loc = children[i]->getAllIndices(indices, loc); + if (is_leaf()) { + for (unsigned int i = 0; i < cum_size; i++) indices[loc + i] = index[i]; + loc += cum_size; + } else { + // Gather indices in children + for (const auto& child : *children) + loc = child.getAllIndices(indices, loc); } return loc; } -template -unsigned int SPTree::getDepth() { - if(is_leaf) return 1; - int depth = 0; - for(unsigned int i = 0; i < no_children; i++) depth = fmax(depth, children[i]->getDepth()); - return 1 + depth; +template +unsigned int SPTreeNode::getDepth() const { + if (is_leaf()) + return 1; + unsigned int depth = 0; + for (const auto& child : *children) depth = max(depth, child.getDepth()); + return 1u + depth; } - // Compute non-edge forces using Barnes-Hut algorithm -template -void SPTree::computeNonEdgeForces(unsigned int point_index, double theta, double neg_f[], double* sum_Q) -{ - +template +void SPTreeNode::computeNonEdgeForces(unsigned int point_index, + const double* data_point, + double neg_f[], double* sum_Q, + double max_width_squared) const { // Make sure that we spend no time on empty nodes or self-interactions - if(cum_size == 0 || (is_leaf && size == 1 && index[0] == point_index)) return; + if (cum_size == 0 || + (cum_size == 1 && __builtin_expect(index[0] == point_index, 0))) + return; // Compute distance between point and center-of-mass + alignas(16) point_t diff; double sqdist = .0; - unsigned int ind = point_index * NDims; - for(unsigned int d = 0; d < NDims; d++) { - buff[d] = data[ind + d] - center_of_mass[d]; - sqdist += buff[d] * buff[d]; + for (int i = 0; i < NDims; ++i) { + diff[i] = data_point[i] - center_of_mass[i]; + sqdist += diff[i] * diff[i]; } // Check whether we can use this node as a "summary" - double max_width = 0.0; - double cur_width; - for(unsigned int d = 0; d < NDims; d++) { - cur_width = boundary.getWidth(d); - max_width = (max_width > cur_width) ? max_width : cur_width; - } - if(is_leaf || max_width / sqrt(sqdist) < theta) { - + if (cum_size == 1) { + sqdist = 1.0 / (1.0 + sqdist); + *sum_Q += sqdist; + double mult = sqdist * sqdist; + for (size_t d = 0; d < NDims; ++d) { + neg_f[d] += mult * diff[d]; + } + } else if (max_width_squared < sqdist) { + // max_width / sqrt(sqdist) < theta // Compute and add t-SNE force between point and current node sqdist = 1.0 / (1.0 + sqdist); double mult = cum_size * sqdist; *sum_Q += mult; mult *= sqdist; - for(unsigned int d = 0; d < NDims; d++) neg_f[d] += mult * buff[d]; + for (size_t d = 0; d < NDims; ++d) { + neg_f[d] += mult * diff[d]; + } + } else { + // Recursively apply Barnes-Hut to children + max_width_squared /= 4.0; + for (const auto& child : *children) + child.computeNonEdgeForces(point_index, data_point, neg_f, sum_Q, + max_width_squared); } - else { +} - // Recursively apply Barnes-Hut to children - for(unsigned int i = 0; i < no_children; i++) children[i]->computeNonEdgeForces(point_index, theta, neg_f, sum_Q); +// Print out tree +template +void SPTreeNode::print(const double* data) const { + if (cum_size == 0) { + fprintf(stderr, "Empty node\n"); + return; + } + + if (is_leaf()) { + fprintf(stderr, "Leaf node; data = ["); + for (int i = 0; i < cum_size; i++) { + const double* point = data + index[i] * NDims; + for (int d = 0; d < NDims; d++) fprintf(stderr, "%f, ", point[d]); + fprintf(stderr, " (index = %d)", index[i]); + if (i < cum_size - 1) + fprintf(stderr, "\n"); + else + fprintf(stderr, "]\n"); + } + } else { + fprintf(stderr, "Intersection node with center-of-mass = ["); + for (const auto& cm : center_of_mass) fprintf(stderr, "%f, ", cm); + fprintf(stderr, "]; children are:\n"); + for (const auto& child : *children) child.print(data); } } +} // namespace _sptree_internal -// Computes edge forces -template -void SPTree::computeEdgeForces(unsigned int* row_P, unsigned int* col_P, double* val_P, int N, double* pos_f) -{ +using namespace _sptree_internal; - // Loop over all edges in the graph - unsigned int ind1 = 0; - unsigned int ind2 = 0; - double sqdist; - for(unsigned int n = 0; n < N; n++) { - for(unsigned int i = row_P[n]; i < row_P[n + 1]; i++) { +template +double SPTree::maxWidth() const { + return *max_element(widths[0].begin(), widths[0].end()); +} + +// Top-node constructor for SPTree -- build tree, too! +template +SPTree::SPTree(const double* inp_data, unsigned int N) : data(inp_data) { + // Compute mean, width, and height of current map (boundaries of SPTree) + int nD = 0; + alignas(16) point_t mean_Y, min_Y, max_Y; + mean_Y.fill(0.0); + min_Y.fill(DBL_MAX); + max_Y.fill(-DBL_MAX); + + const double* elem = inp_data; + for (unsigned int n = 0; n < N; ++n) { + for (int d = 0; d < NDims; ++d) { + mean_Y[d] += elem[d]; + if (elem[d] < min_Y[d]) + min_Y[d] = elem[d]; + if (elem[d] > max_Y[d]) + max_Y[d] = elem[d]; + } + elem += NDims; + } + double dbl_N = static_cast(N); + for (int d = 0; d < NDims; d++) mean_Y[d] /= dbl_N; + + // Construct SPTree + widths.emplace_back(); + point_t& width = widths.back(); + for (int d = 0; d < NDims; d++) { + width[d] = max(max_Y[d] - mean_Y[d], mean_Y[d] - min_Y[d]) + 1e-5; + } + node.setCorner(mean_Y); + fill(N); +} + +// Insert a point into the SPTree +template +bool SPTree::insert(unsigned int new_index) { + return node.insert(new_index, data, &widths, 0); +} + +// Build SPTree on dataset +template +void SPTree::fill(unsigned int N) { + for (unsigned int i = 0; i < N; i++) insert(i); +} + +// Checks whether the specified tree is correct +template +bool SPTree::isCorrect() const { + return node.isCorrect(data, widths.begin()); +} + +// Build a list of all indices in SPTree +template +void SPTree::getAllIndices(unsigned int* indices) const { + node.getAllIndices(indices, 0); +} + +template +unsigned int SPTree::getDepth() const { + return node.getDepth(); +} + +// Compute non-edge forces using Barnes-Hut algorithm +template +void SPTree::computeNonEdgeForces(unsigned int point_index, double theta, + double neg_f[], double* sum_Q) const { + double max_width = maxWidth() / theta; + node.computeNonEdgeForces(point_index, data + point_index * NDims, neg_f, + sum_Q, max_width * max_width); +} + +// Computes edge forces +template +void SPTree::computeEdgeForces(unsigned int* row_P, unsigned int* col_P, + double* val_P, int N, + double* pos_f) const { + // Loop over all edges in the graph + const double* data_1 = data; + for (unsigned int n = 0; n < N; n++) { + for (unsigned int i = row_P[n]; i < row_P[n + 1]; i++) { + double sqdist; // Compute pairwise distance and Q-value sqdist = 1.0; - ind2 = col_P[i] * NDims; + const double* data_2 = data + col_P[i] * NDims; - for(unsigned int d = 0; d < NDims; d++) { - buff[d] = data[ind1 + d] - data[ind2 + d]; - sqdist += buff[d] * buff[d]; + alignas(16) point_t diffs; + for (unsigned int d = 0; d < NDims; d++) { + diffs[d] = data_1[d] - data_2[d]; + sqdist += diffs[d] * diffs[d]; } sqdist = val_P[i] / sqdist; // Sum positive force - for(unsigned int d = 0; d < NDims; d++) pos_f[ind1 + d] += sqdist * buff[d]; + for (unsigned int d = 0; d < NDims; d++) { + pos_f[d] += sqdist * diffs[d]; + } } - ind1 += NDims; + pos_f += NDims; + data_1 += NDims; } } - // Print out tree -template -void SPTree::print() -{ - if(cum_size == 0) { - fprintf(stderr,"Empty node\n"); - return; - } - - if(is_leaf) { - fprintf(stderr,"Leaf node; data = ["); - for(int i = 0; i < size; i++) { - double* point = data + index[i] * NDims; - for(int d = 0; d < NDims; d++) fprintf(stderr,"%f, ", point[d]); - fprintf(stderr," (index = %d)", index[i]); - if(i < size - 1) fprintf(stderr,"\n"); - else fprintf(stderr,"]\n"); - } - } - else { - fprintf(stderr,"Intersection node with center-of-mass = ["); - for(int d = 0; d < NDims; d++) fprintf(stderr,"%f, ", center_of_mass[d]); - fprintf(stderr,"]; children are:\n"); - for(int i = 0; i < no_children; i++) children[i]->print(); - } +template +void SPTree::print() const { + node.print(data); } // declare templates explicitly template class SPTree<2>; template class SPTree<3>; + +} // namespace TSNE diff --git a/tsne/bh_sne_src/sptree.h b/tsne/bh_sne_src/sptree.h index be42686..ba7b325 100644 --- a/tsne/bh_sne_src/sptree.h +++ b/tsne/bh_sne_src/sptree.h @@ -30,93 +30,140 @@ * */ - #ifndef SPTREE_H #define SPTREE_H -using namespace std; +#include +#include +#include + +namespace TSNE { +namespace _sptree_internal { +template +class[[gnu::visibility("internal")]] alignas(16) Cell final { + public: + typedef std::array point_t; -template -class Cell { - double corner[NDims]; - double width[NDims]; + Cell() = default; + explicit Cell(const point_t&); + double getCorner(unsigned int d) const; + void setCorner(const point_t& inp_corner); + void setCorner(unsigned int d, double val); + bool containsPoint(const double* point, const point_t& width) const; -public: - Cell(); - Cell(double* inp_corner, double* inp_width); - ~Cell(); + private: + point_t corner; - double getCorner(unsigned int d); - double getWidth(unsigned int d); - void setCorner(unsigned int d, double val); - void setWidth(unsigned int d, double val); - bool containsPoint(double point[]); + // disallow copy + Cell(const Cell&) = delete; }; -template -class SPTree -{ -public: - enum { no_children = 2 * SPTree::no_children }; - -private: - // Fixed constants - static const unsigned int QT_NODE_CAPACITY = 1; - - // A buffer we use when doing force computations - double buff[NDims]; - - // Properties of this node in the tree - SPTree* parent; - unsigned int dimension; - bool is_leaf; - unsigned int size; - unsigned int cum_size; - - // Axis-aligned bounding box stored as a center with half-dimensions to represent the boundaries of this quad tree - Cell boundary; - - // Indices in this space-partitioning tree node, corresponding center-of-mass, and list of all children - double* data; - double center_of_mass[NDims]; - unsigned int index[QT_NODE_CAPACITY]; - - // Children - SPTree* children[no_children]; - -public: - SPTree(double* inp_data, unsigned int N); - SPTree(double* inp_data, double* inp_corner, double* inp_width); - SPTree(double* inp_data, unsigned int N, double* inp_corner, double* inp_width); - SPTree(SPTree* inp_parent, double* inp_data, unsigned int N, double* inp_corner, double* inp_width); - SPTree(SPTree* inp_parent, double* inp_data, double* inp_corner, double* inp_width); - ~SPTree(); - void setData(double* inp_data); - SPTree* getParent(); - void construct(Cell boundary); - bool insert(unsigned int new_index); - void subdivide(); - bool isCorrect(); - void rebuildTree(); - void getAllIndices(unsigned int* indices); - unsigned int getDepth(); - void computeNonEdgeForces(unsigned int point_index, double theta, double neg_f[], double* sum_Q); - void computeEdgeForces(unsigned int* row_P, unsigned int* col_P, double* val_P, int N, double* pos_f); - void print(); - -private: - void init(SPTree* inp_parent, double* inp_data, double* inp_corner, double* inp_width); - void fill(unsigned int N); - unsigned int getAllIndices(unsigned int* indices, unsigned int loc); - bool isChild(unsigned int test_index, unsigned int start, unsigned int end); +template +class[[gnu::visibility("internal")]] alignas(16) SPTreeNode final { + public: + typedef typename Cell::point_t point_t; + enum { no_children = 2 * SPTreeNode::no_children }; + + private: + // Fixed constants + static constexpr unsigned int QT_NODE_CAPACITY = 1; + + // Axis-aligned bounding box stored as a center with half-dimensions to + // represent the boundaries of this quad tree + point_t center_of_mass; + Cell boundary; + + // Children + std::unique_ptr, no_children>> children; + + // Properties of this node in the tree + unsigned int cum_size; + + // Indices in this space-partitioning tree node, corresponding + // center-of-mass, and list of all children + std::array index; + + // Disallow copy + SPTreeNode(const SPTreeNode&) = delete; + + void subdivide(const double* data, std::vector* widths, + typename std::vector::size_type depth); + unsigned int which_child(const double* point) const; + void make_child(unsigned int i, const point_t& width); + + bool is_leaf() const; + + public: + SPTreeNode(); + explicit SPTreeNode(const point_t& corner); + + void setCorner(const point_t& corner); + + bool insert(unsigned int new_index, const double* data, + std::vector* widths, + typename std::vector::size_type depth); + void computeNonEdgeForces[[gnu::hot]]( + unsigned int point_index, const double* data_point, double neg_f[], + double* sum_Q, double max_width_squared) const; + + // Methods used for debugging + bool isCorrect[[gnu::cold]]( + const double* data, typename std::vector::const_iterator width) + const; + + unsigned int getAllIndices[[gnu::cold]](unsigned int* indices, + unsigned int loc) const; + unsigned int getDepth[[gnu::cold]]() const; + void print[[gnu::cold]](const double* data) const; }; template <> -struct SPTree<0> -{ - enum { no_children = 1 }; +class SPTreeNode<0> { + public: + enum { no_children = 1 }; }; +} // namespace _sptree_internal + +template +class [[gnu::visibility("internal")]] SPTree { + public: + typedef typename _sptree_internal::SPTreeNode::point_t point_t; + enum { no_children = _sptree_internal::SPTreeNode::no_children }; + + private: + _sptree_internal::SPTreeNode node; + + // The width for each cell is the same at each level. The top node owns + // this and the children get references to it. + std::vector widths; + bool insert(unsigned int new_index); + + const double* data; + double maxWidth() const; + + public: + SPTree(const double* inp_data, unsigned int N); + + void computeNonEdgeForces[[gnu::hot]](unsigned int point_index, double theta, + double neg_f[], double* sum_Q) const; + void computeEdgeForces[[gnu::hot]](unsigned int* row_P, unsigned int* col_P, + double* val_P, int N, double* pos_f) const; + + // methods used for debugging. + bool isCorrect[[gnu::cold]]() const; + void getAllIndices[[gnu::cold]](unsigned int* indices) const; + unsigned int getDepth[[gnu::cold]]() const; + void print[[gnu::cold]]() const; + + private: + void fill(unsigned int N); + // Disallow copy + SPTree(const SPTree&) = delete; +}; + +} // namespace TSNE + #endif diff --git a/tsne/bh_sne_src/tsne.cpp b/tsne/bh_sne_src/tsne.cpp index 5fac463..8f2368e 100644 --- a/tsne/bh_sne_src/tsne.cpp +++ b/tsne/bh_sne_src/tsne.cpp @@ -30,51 +30,91 @@ * */ -#include +#include "tsne.h" + +#include +#include +#include +#include +#include +#include +#include #include -#include +#include #include -#include -#include -#include -#include -#include -#include -#include "vptree.h" + #include "sptree.h" -#include "tsne.h" +#include "vptree.h" -using namespace std; +using std::array; +using std::move; +using std::vector; -#ifdef TSNE3D -#define NDIMS 3 -#else -#define NDIMS 2 -#endif +namespace TSNE { +namespace { -// Perform t-SNE -void TSNE::run(double* X, int N, int D, double* Y, int no_dims, double perplexity, double theta, int rand_seed, - bool skip_random_init, double *init, bool use_init, - int max_iter, int stop_lying_iter, int mom_switch_iter - ) { +static inline double sign(double x) { + return (x == .0 ? .0 : (x < .0 ? -1.0 : 1.0)); +} + +void symmetrizeMatrix(vector* row_P, vector* col_P, + vector* val_P, int N); +template +void computeGradient(double* P, unsigned int* inp_row_P, + unsigned int* inp_col_P, double* inp_val_P, + const double* Y, int N, double* dC, double theta); +template +void computeExactGradient(double* P, const double* Y, int N, double* dC); +template +double evaluateError(double* P, const double* Y, int N); +template +double evaluateError(unsigned int* row_P, unsigned int* col_P, double* val_P, + const double* Y, int N, double theta); +void zeroMean(double* X, int N, int D); +template +void zeroMean(double* X, int N); +void computeGaussianPerplexity(const double* X, int N, int D, double* P, + double perplexity); +void computeGaussianPerplexity(const double* X, int N, int D, + vector* _row_P, + vector* _col_P, + vector* _val_P, double perplexity, + int K); +template +void computeSquaredEuclideanDistance(const double* X, int N, double* DD); +void computeSquaredEuclideanDistance(const double* X, int N, int D, double* DD); +double randn(); +// Perform t-SNE +template +void run(double* X, int N, int D, double* Y, double perplexity, double theta, + int rand_seed, bool skip_random_init, double* init, bool use_init, + int max_iter, int stop_lying_iter, int mom_switch_iter) { // Set random seed if (skip_random_init != true) { - if(rand_seed >= 0) { - fprintf(stderr,"Using random seed: %d\n", rand_seed); - srand((unsigned int) rand_seed); - } else { - fprintf(stderr,"Using current time as random seed...\n"); - srand(time(NULL)); - } + if (rand_seed >= 0) { + fprintf(stderr, "Using random seed: %d\n", rand_seed); + srand((unsigned int)rand_seed); + } else { + fprintf(stderr, "Using current time as random seed...\n"); + srand(time(NULL)); + } } // Determine whether we are using an exact algorithm - if(N - 1 < 3 * perplexity) { fprintf(stderr,"Perplexity too large for the number of data points!\n"); exit(1); } - fprintf(stderr,"Using no_dims = %d, perplexity = %f, and theta = %f\n", no_dims, perplexity, theta); - bool exact = (theta == .0) ? true : false; + if (N - 1 < 3 * perplexity) { + fprintf(stderr, + "Perplexity too large for the number of data points!\n"); + exit(1); + } + fprintf(stderr, + "Using D = %d, no_dims = %d, perplexity = %f, and theta = %f\n", D, + no_dims, perplexity, theta); + const bool exact = __builtin_expect((theta == .0) ? true : false, 0); - fprintf(stderr,"Using max_iter = %d, stop_lying_iter = %d, mom_switch_iter = %d\n", max_iter, stop_lying_iter, mom_switch_iter); + fprintf(stderr, + "Using max_iter = %d, stop_lying_iter = %d, mom_switch_iter = %d\n", + max_iter, stop_lying_iter, mom_switch_iter); // Set learning parameters float total_time = .0; @@ -83,182 +123,198 @@ void TSNE::run(double* X, int N, int D, double* Y, int no_dims, double perplexit double eta = 200.0; // Allocate some memory - double* dY = (double*) malloc(N * no_dims * sizeof(double)); - double* uY = (double*) malloc(N * no_dims * sizeof(double)); - double* gains = (double*) malloc(N * no_dims * sizeof(double)); - if(dY == NULL || uY == NULL || gains == NULL) { fprintf(stderr,"Memory allocation failed!\n"); exit(1); } - for(int i = 0; i < N * no_dims; i++) uY[i] = .0; - for(int i = 0; i < N * no_dims; i++) gains[i] = 1.0; + vector dY(N * no_dims); + vector uY(N * no_dims, .0); + vector gains(N * no_dims, 1.0); // Normalize input data (to prevent numerical problems) - fprintf(stderr,"Computing input similarities...\n"); + fprintf(stderr, "Computing input similarities...\n"); start = clock(); zeroMean(X, N, D); double max_X = .0; - for(int i = 0; i < N * D; i++) { - if(fabs(X[i]) > max_X) max_X = fabs(X[i]); + for (int i = 0; i < N * D; i++) { + if (fabs(X[i]) > max_X) + max_X = fabs(X[i]); } - for(int i = 0; i < N * D; i++) X[i] /= max_X; + for (int i = 0; i < N * D; i++) X[i] /= max_X; // Compute input similarities for exact t-SNE - double* P; unsigned int* row_P; unsigned int* col_P; double* val_P; - if(exact) { - + vector P; + vector row_P, col_P; + vector val_P; + if (exact) { // Compute similarities - fprintf(stderr,"Exact?"); - P = (double*) malloc(N * N * sizeof(double)); - if(P == NULL) { fprintf(stderr,"Memory allocation failed!\n"); exit(1); } - computeGaussianPerplexity(X, N, D, P, perplexity); + P.resize(N * N); + fprintf(stderr, "Computing exact perplexity...\n"); + computeGaussianPerplexity(X, N, D, P.data(), perplexity); // Symmetrize input similarities - fprintf(stderr,"Symmetrizing...\n"); + fprintf(stderr, "Symmetrizing...\n"); int nN = 0; - for(int n = 0; n < N; n++) { + for (int n = 0; n < N; n++) { int mN = (n + 1) * N; - for(int m = n + 1; m < N; m++) { + for (int m = n + 1; m < N; m++) { P[nN + m] += P[mN + n]; - P[mN + n] = P[nN + m]; + P[mN + n] = P[nN + m]; mN += N; } nN += N; } double sum_P = .0; - for(int i = 0; i < N * N; i++) sum_P += P[i]; - for(int i = 0; i < N * N; i++) P[i] /= sum_P; + for (int i = 0; i < N * N; i++) sum_P += P[i]; + for (int i = 0; i < N * N; i++) P[i] /= sum_P; } // Compute input similarities for approximate t-SNE else { - + fprintf(stderr, "Computing approximate perplexity...\n"); // Compute asymmetric pairwise input similarities - computeGaussianPerplexity(X, N, D, &row_P, &col_P, &val_P, perplexity, (int) (3 * perplexity)); + computeGaussianPerplexity(X, N, D, &row_P, &col_P, &val_P, perplexity, + (int)(3 * perplexity)); // Symmetrize input similarities symmetrizeMatrix(&row_P, &col_P, &val_P, N); double sum_P = .0; - for(int i = 0; i < row_P[N]; i++) sum_P += val_P[i]; - for(int i = 0; i < row_P[N]; i++) val_P[i] /= sum_P; + for (int i = 0; i < row_P[N]; i++) sum_P += val_P[i]; + for (int i = 0; i < row_P[N]; i++) val_P[i] /= sum_P; } end = clock(); // Lie about the P-values - if(exact) { for(int i = 0; i < N * N; i++) P[i] *= 12.0; } - else { for(int i = 0; i < row_P[N]; i++) val_P[i] *= 12.0; } + if (exact) { + for (int i = 0; i < N * N; i++) P[i] *= 12.0; + } else { + for (int i = 0; i < row_P[N]; i++) val_P[i] *= 12.0; + } // Initialize solution (randomly or with given coordinates) if (!use_init && !skip_random_init) { - for(int i = 0; i < N * no_dims; i++) { - Y[i] = randn() * .0001; - } + for (int i = 0; i < N * no_dims; i++) { + Y[i] = randn() * .0001; + } } else if (use_init) { - for(int i = 0; i < N * no_dims; i++) { - Y[i] = init[i]; - } + for (int i = 0; i < N * no_dims; i++) { + Y[i] = init[i]; + } } - // Perform main training loop - if(exact) fprintf(stderr,"Input similarities computed in %4.2f seconds!\nLearning embedding...\n", (float) (end - start) / CLOCKS_PER_SEC); - else fprintf(stderr,"Input similarities computed in %4.2f seconds (sparsity = %f)!\nLearning embedding...\n", (float) (end - start) / CLOCKS_PER_SEC, (double) row_P[N] / ((double) N * (double) N)); + // Perform main training loop + if (exact) + fprintf(stderr, + "Input similarities computed in %4.2f seconds!\nLearning " + "embedding...\n", + (float)(end - start) / CLOCKS_PER_SEC); + else + fprintf(stderr, + "Input similarities computed in %4.2f seconds (sparsity = " + "%f)!\nLearning embedding...\n", + (float)(end - start) / CLOCKS_PER_SEC, + (double)row_P[N] / ((double)N * (double)N)); start = clock(); - for(int iter = 0; iter < max_iter; iter++) { - + for (int iter = 0; iter < max_iter; iter++) { // Compute (approximate) gradient - if(exact) computeExactGradient(P, Y, N, no_dims, dY); - else computeGradient(P, row_P, col_P, val_P, Y, N, no_dims, dY, theta); + if (exact) + computeExactGradient(P.data(), Y, N, dY.data()); + else + computeGradient(P.data(), row_P.data(), col_P.data(), + val_P.data(), Y, N, dY.data(), theta); // Update gains - for(int i = 0; i < N * no_dims; i++) gains[i] = (sign(dY[i]) != sign(uY[i])) ? (gains[i] + .2) : (gains[i] * .8); - for(int i = 0; i < N * no_dims; i++) if(gains[i] < .01) gains[i] = .01; + for (int i = 0; i < N * no_dims; i++) + gains[i] = (sign(dY[i]) != sign(uY[i])) ? (gains[i] + .2) + : (gains[i] * .8); + for (int i = 0; i < N * no_dims; i++) + if (gains[i] < .01) + gains[i] = .01; // Perform gradient update (with momentum and gains) - for(int i = 0; i < N * no_dims; i++) uY[i] = momentum * uY[i] - eta * gains[i] * dY[i]; - for(int i = 0; i < N * no_dims; i++) Y[i] = Y[i] + uY[i]; + for (int i = 0; i < N * no_dims; i++) + uY[i] = momentum * uY[i] - eta * gains[i] * dY[i]; + for (int i = 0; i < N * no_dims; i++) Y[i] = Y[i] + uY[i]; // Make solution zero-mean - zeroMean(Y, N, no_dims); + zeroMean(Y, N); // Stop lying about the P-values after a while, and switch momentum - if(iter == stop_lying_iter) { - if(exact) { for(int i = 0; i < N * N; i++) P[i] /= 12.0; } - else { for(int i = 0; i < row_P[N]; i++) val_P[i] /= 12.0; } + if (iter == stop_lying_iter) { + if (exact) { + for (int i = 0; i < N * N; i++) P[i] /= 12.0; + } else { + for (int i = 0; i < row_P[N]; i++) val_P[i] /= 12.0; + } } - if(iter == mom_switch_iter) momentum = final_momentum; + if (iter == mom_switch_iter) + momentum = final_momentum; // Print out progress if (iter > 0 && (iter % 50 == 0 || iter == max_iter - 1)) { end = clock(); double C = .0; - if(exact) C = evaluateError(P, Y, N, no_dims); - else C = evaluateError(row_P, col_P, val_P, Y, N, no_dims, theta); // doing approximate computation here! - if(iter == 0) - fprintf(stderr,"Iteration %d: error is %f\n", iter + 1, C); + if (exact) + C = evaluateError(P.data(), Y, N); + else + C = evaluateError( + row_P.data(), col_P.data(), val_P.data(), Y, N, + theta); // doing approximate computation here! + if (iter == 0) + fprintf(stderr, "Iteration %d: error is %f\n", iter + 1, C); else { - total_time += (float) (end - start) / CLOCKS_PER_SEC; - fprintf(stderr,"Iteration %d: error is %f (50 iterations in %4.2f seconds)\n", iter, C, (float) (end - start) / CLOCKS_PER_SEC); + total_time += (float)(end - start) / CLOCKS_PER_SEC; + fprintf(stderr, + "Iteration %d: error is %f (50 iterations in %4.2f " + "seconds)\n", + iter, C, (float)(end - start) / CLOCKS_PER_SEC); } - start = clock(); + start = clock(); } } - end = clock(); total_time += (float) (end - start) / CLOCKS_PER_SEC; + end = clock(); + total_time += (float)(end - start) / CLOCKS_PER_SEC; - // Clean up memory - free(dY); - free(uY); - free(gains); - if(exact) free(P); - else { - free(row_P); row_P = NULL; - free(col_P); col_P = NULL; - free(val_P); val_P = NULL; - } - fprintf(stderr,"Fitting performed in %4.2f seconds.\n", total_time); + fprintf(stderr, "Fitting performed in %4.2f seconds.\n", total_time); } - // Compute gradient of the t-SNE cost function (using Barnes-Hut algorithm) -void TSNE::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) -{ - +template +void computeGradient(double* P, unsigned int* inp_row_P, + unsigned int* inp_col_P, double* inp_val_P, + const double* Y, int N, double* dC, double theta) { // Construct space-partitioning tree on current map - SPTree* tree = new SPTree(Y, N); + SPTree tree(Y, N); // Compute all terms required for t-SNE gradient double sum_Q = .0; - double* pos_f = (double*) calloc(N * D, sizeof(double)); - double* neg_f = (double*) calloc(N * D, sizeof(double)); - if(pos_f == NULL || neg_f == NULL) { fprintf(stderr,"Memory allocation failed!\n"); exit(1); } - tree->computeEdgeForces(inp_row_P, inp_col_P, inp_val_P, N, pos_f); - for(int n = 0; n < N; n++) tree->computeNonEdgeForces(n, theta, neg_f + n * D, &sum_Q); + auto len = N * D; + vector pos_f(2 * len); + double* neg_f = pos_f.data() + len; + tree.computeEdgeForces(inp_row_P, inp_col_P, inp_val_P, N, pos_f.data()); + for (int n = 0; n < N; n++) + tree.computeNonEdgeForces(n, theta, neg_f + n * D, &sum_Q); // Compute final t-SNE gradient - for(int i = 0; i < N * D; i++) { + for (int i = 0; i < len; i++) { dC[i] = pos_f[i] - (neg_f[i] / sum_Q); } - free(pos_f); - free(neg_f); - delete tree; } // Compute gradient of the t-SNE cost function (exact) -void TSNE::computeExactGradient(double* P, double* Y, int N, int D, double* dC) { - - // Make sure the current gradient contains zeros - for(int i = 0; i < N * D; i++) dC[i] = 0.0; +template +void computeExactGradient(double* P, const double* Y, int N, double* dC) { + // Make sure the current gradient contains zeros + for (int i = 0; i < N * D; i++) dC[i] = 0.0; // Compute the squared Euclidean distance matrix - double* DD = (double*) malloc(N * N * sizeof(double)); - if(DD == NULL) { fprintf(stderr,"Memory allocation failed!\n"); exit(1); } - computeSquaredEuclideanDistance(Y, N, D, DD); + vector DD(N * N); + computeSquaredEuclideanDistance(Y, N, DD.data()); // Compute Q-matrix and normalization sum - double* Q = (double*) malloc(N * N * sizeof(double)); - if(Q == NULL) { fprintf(stderr,"Memory allocation failed!\n"); exit(1); } + vector Q(N * N); double sum_Q = .0; int nN = 0; - for(int n = 0; n < N; n++) { - for(int m = 0; m < N; m++) { - if(n != m) { + 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]); sum_Q += Q[nN + m]; } @@ -266,244 +322,228 @@ void TSNE::computeExactGradient(double* P, double* Y, int N, int D, double* dC) nN += N; } - // Perform the computation of the gradient + // Perform the computation of the gradient nN = 0; int nD = 0; - for(int n = 0; n < N; n++) { + for (int n = 0; n < N; n++) { int mD = 0; - for(int m = 0; m < N; m++) { - if(n != m) { + for (int m = 0; m < N; m++) { + if (n != m) { double mult = (P[nN + m] - (Q[nN + m] / sum_Q)) * Q[nN + m]; - for(int d = 0; d < D; d++) { + for (int d = 0; d < D; d++) { dC[nD + d] += (Y[nD + d] - Y[mD + d]) * mult; } } mD += D; - } + } nN += N; nD += D; - } - - // Free memory - free(DD); DD = NULL; - free(Q); Q = NULL; + } } - // Evaluate t-SNE cost function (exactly) -double TSNE::evaluateError(double* P, double* Y, int N, int D) { - +template +double evaluateError(double* P, const double* Y, int N) { // Compute the squared Euclidean distance matrix - double* DD = (double*) malloc(N * N * sizeof(double)); - double* Q = (double*) malloc(N * N * sizeof(double)); - if(DD == NULL || Q == NULL) { fprintf(stderr,"Memory allocation failed!\n"); exit(1); } - computeSquaredEuclideanDistance(Y, N, D, DD); + vector DD(N * N); + computeSquaredEuclideanDistance(Y, N, DD.data()); // Compute Q-matrix and normalization sum int nN = 0; double sum_Q = DBL_MIN; - for(int n = 0; n < N; n++) { - for(int m = 0; m < N; m++) { - if(n != m) { + vector Q(N * N); + 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]); sum_Q += Q[nN + m]; - } - else Q[nN + m] = DBL_MIN; + } else + Q[nN + m] = DBL_MIN; } nN += N; } - for(int i = 0; i < N * N; i++) Q[i] /= sum_Q; + for (int i = 0; i < N * N; i++) Q[i] /= sum_Q; // Sum t-SNE error double C = .0; - for(int n = 0; n < N * N; n++) { - C += P[n] * log((P[n] + FLT_MIN) / (Q[n] + FLT_MIN)); + for (int n = 0; n < N * N; n++) { + C += P[n] * log((P[n] + FLT_MIN) / (Q[n] + FLT_MIN)); } - // Clean up memory - free(DD); - free(Q); - return C; + return C; } // Evaluate t-SNE cost function (approximately) -double TSNE::evaluateError(unsigned int* row_P, unsigned int* col_P, double* val_P, double* Y, int N, int D, double theta) -{ - +template +double evaluateError(unsigned int* row_P, unsigned int* col_P, double* val_P, + const double* Y, int N, double theta) { // Get estimate of normalization term - SPTree* tree = new SPTree(Y, N); - double* buff = (double*) calloc(D, sizeof(double)); + SPTree tree(Y, N); + double buff[D]; double sum_Q = .0; - for(int n = 0; n < N; n++) tree->computeNonEdgeForces(n, theta, buff, &sum_Q); + for (int n = 0; n < N; n++) + tree.computeNonEdgeForces(n, theta, buff, &sum_Q); // Loop over all edges to compute t-SNE error int ind1, ind2; double C = .0, Q; - for(int n = 0; n < N; n++) { + for (int n = 0; n < N; n++) { ind1 = n * D; - for(int i = row_P[n]; i < row_P[n + 1]; i++) { + for (int i = row_P[n]; i < row_P[n + 1]; i++) { Q = .0; ind2 = col_P[i] * D; - 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]; + 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; C += val_P[i] * log((val_P[i] + FLT_MIN) / (Q + FLT_MIN)); } } - - // Clean up memory - free(buff); - delete tree; return C; } - // Compute input similarities with a fixed perplexity -void TSNE::computeGaussianPerplexity(double* X, int N, int D, double* P, double perplexity) { - - // Compute the squared Euclidean distance matrix - double* DD = (double*) malloc(N * N * sizeof(double)); - if(DD == NULL) { fprintf(stderr,"Memory allocation failed!\n"); exit(1); } - computeSquaredEuclideanDistance(X, N, D, DD); +void computeGaussianPerplexity(const double* X, int N, int D, double* P, + double perplexity) { + // Compute the squared Euclidean distance matrix + vector DD(N * N); + computeSquaredEuclideanDistance(X, N, D, DD.data()); - // Compute the Gaussian kernel row by row - int nN = 0; - for(int n = 0; n < N; n++) { - - // Initialize some variables - bool found = false; - double beta = 1.0; - double min_beta = -DBL_MAX; - double max_beta = DBL_MAX; - double tol = 1e-5; + // Compute the Gaussian kernel row by row + double* rP = P; + double* rD = DD.data(); + for (int n = 0; n < N; n++) { + // Initialize some variables + bool found = false; + double beta = 1.0; + double min_beta = -DBL_MAX; + double max_beta = DBL_MAX; + double tol = 1e-5; double sum_P; - // Iterate until we found a good perplexity - int iter = 0; - while(!found && iter < 200) { - - // Compute Gaussian kernel row - for(int m = 0; m < N; m++) P[nN + m] = exp(-beta * DD[nN + m]); - P[nN + n] = DBL_MIN; - - // Compute entropy of current row - sum_P = DBL_MIN; - for(int m = 0; m < N; m++) sum_P += P[nN + m]; - double H = 0.0; - for(int m = 0; m < N; m++) H += beta * (DD[nN + m] * P[nN + m]); - H = (H / sum_P) + log(sum_P); - - // Evaluate whether the entropy is within the tolerance level - double Hdiff = H - log(perplexity); - if(Hdiff < tol && -Hdiff < tol) { - found = true; - } - else { - if(Hdiff > 0) { - min_beta = beta; - if(max_beta == DBL_MAX || max_beta == -DBL_MAX) - beta *= 2.0; - else - beta = (beta + max_beta) / 2.0; - } - else { - max_beta = beta; - if(min_beta == -DBL_MAX || min_beta == DBL_MAX) - beta /= 2.0; - else - beta = (beta + min_beta) / 2.0; - } - } - - // Update iteration counter - iter++; - } - - // Row normalize P - for(int m = 0; m < N; m++) P[nN + m] /= sum_P; - nN += N; - } + // Iterate until we found a good perplexity + int iter = 0; + while (!found && iter < 200) { + // Compute Gaussian kernel row + for (int m = 0; m < N; m++) rP[m] = exp(-beta * rD[m]); + rP[n] = DBL_MIN; - // Clean up memory - free(DD); DD = NULL; -} + // Compute entropy of current row + sum_P = DBL_MIN; + for (int m = 0; m < N; m++) sum_P += rP[m]; + double H = 0.0; + for (int m = 0; m < N; m++) H += rD[m] * rP[m]; + H *= beta; + H = (H / sum_P) + log(sum_P); + + // Evaluate whether the entropy is within the tolerance level + double Hdiff = H - log(perplexity); + if (Hdiff < tol && -Hdiff < tol) { + found = true; + } else { + if (Hdiff > 0) { + min_beta = beta; + if (max_beta == DBL_MAX || max_beta == -DBL_MAX) + beta *= 2.0; + else + beta = (beta + max_beta) / 2.0; + } else { + max_beta = beta; + if (min_beta == -DBL_MAX || min_beta == DBL_MAX) + beta /= 2.0; + else + beta = (beta + min_beta) / 2.0; + } + } + // Update iteration counter + iter++; + } -// Compute input similarities with a fixed perplexity using ball trees (this function allocates memory another function should free) -void TSNE::computeGaussianPerplexity(double* X, int N, int D, unsigned int** _row_P, unsigned int** _col_P, double** _val_P, double perplexity, int K) { + // Row normalize P + for (int m = 0; m < N; m++) rP[m] /= sum_P; + rP += N; + rD += N; + } +} - if(perplexity > K) fprintf(stderr,"Perplexity should be lower than K!\n"); +// Compute input similarities with a fixed perplexity using ball trees (this +// function allocates memory another function should free) +void computeGaussianPerplexity(const double* X, int N, int D, + vector* _row_P, + vector* _col_P, + vector* _val_P, double perplexity, + int K) { + if (perplexity > K) + fprintf(stderr, "Perplexity should be lower than K!\n"); // Allocate the memory we need - *_row_P = (unsigned int*) malloc((N + 1) * sizeof(unsigned int)); - *_col_P = (unsigned int*) calloc(N * K, sizeof(unsigned int)); - *_val_P = (double*) calloc(N * K, sizeof(double)); - if(*_row_P == NULL || *_col_P == NULL || *_val_P == NULL) { fprintf(stderr,"Memory allocation failed!\n"); exit(1); } - unsigned int* row_P = *_row_P; - unsigned int* col_P = *_col_P; - double* val_P = *_val_P; - double* cur_P = (double*) malloc((N - 1) * sizeof(double)); - if(cur_P == NULL) { fprintf(stderr,"Memory allocation failed!\n"); exit(1); } + _row_P->resize(N + 1); + _col_P->resize(N * K, 0); + _val_P->resize(N * K, 0); + vector& row_P = *_row_P; + vector& col_P = *_col_P; + vector& val_P = *_val_P; + vector cur_P(N - 1); row_P[0] = 0; - for(int n = 0; n < N; n++) row_P[n + 1] = row_P[n] + (unsigned int) K; + for (int n = 0; n < N; n++) row_P[n + 1] = row_P[n] + (unsigned int)K; // Build ball tree on data set - VpTree* tree = new VpTree(); - vector obj_X(N, DataPoint(D, -1, X)); - for(int n = 0; n < N; n++) obj_X[n] = DataPoint(D, n, X + n * D); - tree->create(obj_X); + VpTree tree((euclidean_distance(D))); + vector obj_X(N); + for (int n = 0; n < N; n++) obj_X[n] = DataPoint(X + n * D); + tree.create(obj_X); // Loop over all points to find nearest neighbors - fprintf(stderr,"Building tree...\n"); + fprintf(stderr, "Building tree...\n"); vector indices; vector distances; - for(int n = 0; n < N; n++) { - - if(n % 10000 == 0) fprintf(stderr," - point %d of %d\n", n, N); + for (int n = 0; n < N; n++) { + if (n % 10000 == 0) + fprintf(stderr, " - point %d of %d\n", n, N); // Find nearest neighbors indices.clear(); distances.clear(); - tree->search(obj_X[n], K + 1, &indices, &distances); + tree.search(obj_X[n], K + 1, &indices, &distances); // Initialize some variables for binary search bool found = false; double beta = 1.0; double min_beta = -DBL_MAX; - double max_beta = DBL_MAX; + double max_beta = DBL_MAX; double tol = 1e-5; // Iterate until we found a good perplexity - int iter = 0; double sum_P; - while(!found && iter < 200) { - + int iter = 0; + double sum_P; + while (!found && iter < 200) { // Compute Gaussian kernel row - for(int m = 0; m < K; m++) cur_P[m] = exp(-beta * distances[m + 1] * distances[m + 1]); + for (int m = 0; m < K; m++) + cur_P[m] = exp(-beta * distances[m + 1] * distances[m + 1]); // Compute entropy of current row sum_P = DBL_MIN; - for(int m = 0; m < K; m++) sum_P += cur_P[m]; + for (int m = 0; m < K; m++) sum_P += cur_P[m]; double H = .0; - for(int m = 0; m < K; m++) H += beta * (distances[m + 1] * distances[m + 1] * cur_P[m]); - H = (H / sum_P) + log(sum_P); + for (int m = 0; m < K; m++) + H += (distances[m + 1] * distances[m + 1] * cur_P[m]); + H = (H * beta / sum_P) + log(sum_P); // Evaluate whether the entropy is within the tolerance level double Hdiff = H - log(perplexity); - if(Hdiff < tol && -Hdiff < tol) { + if (Hdiff < tol && -Hdiff < tol) { found = true; - } - else { - if(Hdiff > 0) { + } else { + if (Hdiff > 0) { min_beta = beta; - if(max_beta == DBL_MAX || max_beta == -DBL_MAX) + if (max_beta == DBL_MAX || max_beta == -DBL_MAX) beta *= 2.0; else beta = (beta + max_beta) / 2.0; - } - else { + } else { max_beta = beta; - if(min_beta == -DBL_MAX || min_beta == DBL_MAX) + if (min_beta == -DBL_MAX || min_beta == DBL_MAX) beta /= 2.0; else beta = (beta + min_beta) / 2.0; @@ -515,40 +555,35 @@ void TSNE::computeGaussianPerplexity(double* X, int N, int D, unsigned int** _ro } // Row-normalize current row of P and store in matrix - for(unsigned int m = 0; m < K; m++) cur_P[m] /= sum_P; - for(unsigned int m = 0; m < K; m++) { - col_P[row_P[n] + m] = (unsigned int) indices[m + 1].index(); + for (unsigned int m = 0; m < K; m++) cur_P[m] /= sum_P; + for (unsigned int m = 0; m < K; m++) { + col_P[row_P[n] + m] = (unsigned int)(indices[m + 1]._x - X) / D; val_P[row_P[n] + m] = cur_P[m]; } } - - // Clean up memory - obj_X.clear(); - free(cur_P); - delete tree; } - // Symmetrizes a sparse matrix -void TSNE::symmetrizeMatrix(unsigned int** _row_P, unsigned int** _col_P, double** _val_P, int N) { - +void symmetrizeMatrix(vector* _row_P, + vector* _col_P, vector* _val_P, + int N) { // Get sparse matrix - unsigned int* row_P = *_row_P; - unsigned int* col_P = *_col_P; - double* val_P = *_val_P; + vector& row_P = *_row_P; + vector& col_P = *_col_P; + vector& val_P = *_val_P; // Count number of elements and row counts of symmetric matrix - int* row_counts = (int*) calloc(N, sizeof(int)); - if(row_counts == NULL) { fprintf(stderr,"Memory allocation failed!\n"); exit(1); } - for(int n = 0; n < N; n++) { - for(int i = row_P[n]; i < row_P[n + 1]; i++) { - + vector row_counts(N); + for (int n = 0; n < N; n++) { + for (int i = row_P[n]; i < row_P[n + 1]; i++) { // Check whether element (col_P[i], n) is present bool present = false; - for(int m = row_P[col_P[i]]; m < row_P[col_P[i] + 1]; m++) { - if(col_P[m] == n) present = true; + for (int m = row_P[col_P[i]]; m < row_P[col_P[i] + 1]; m++) { + if (col_P[m] == n) + present = true; } - if(present) row_counts[n]++; + if (present) + row_counts[n]++; else { row_counts[n]++; row_counts[col_P[i]]++; @@ -556,78 +591,80 @@ void TSNE::symmetrizeMatrix(unsigned int** _row_P, unsigned int** _col_P, double } } int no_elem = 0; - for(int n = 0; n < N; n++) no_elem += row_counts[n]; + for (int n = 0; n < N; n++) no_elem += row_counts[n]; // Allocate memory for symmetrized matrix - unsigned int* sym_row_P = (unsigned int*) malloc((N + 1) * sizeof(unsigned int)); - unsigned int* sym_col_P = (unsigned int*) malloc(no_elem * sizeof(unsigned int)); - double* sym_val_P = (double*) malloc(no_elem * sizeof(double)); - if(sym_row_P == NULL || sym_col_P == NULL || sym_val_P == NULL) { fprintf(stderr,"Memory allocation failed!\n"); exit(1); } + vector sym_row_P(N + 1); + vector sym_col_P(no_elem); + vector sym_val_P(no_elem); // Construct new row indices for symmetric matrix sym_row_P[0] = 0; - for(int n = 0; n < N; n++) sym_row_P[n + 1] = sym_row_P[n] + (unsigned int) row_counts[n]; + for (int n = 0; n < N; n++) + sym_row_P[n + 1] = sym_row_P[n] + (unsigned int)row_counts[n]; // Fill the result matrix - int* offset = (int*) calloc(N, sizeof(int)); - if(offset == NULL) { fprintf(stderr,"Memory allocation failed!\n"); exit(1); } - for(int n = 0; n < N; n++) { - for(unsigned int i = row_P[n]; i < row_P[n + 1]; i++) { // considering element(n, col_P[i]) + vector offset(N); + for (int n = 0; n < N; n++) { + for (unsigned int i = row_P[n]; i < row_P[n + 1]; + i++) { // considering element(n, col_P[i]) // Check whether element (col_P[i], n) is present bool present = false; - for(unsigned int m = row_P[col_P[i]]; m < row_P[col_P[i] + 1]; m++) { - if(col_P[m] == n) { + for (unsigned int m = row_P[col_P[i]]; m < row_P[col_P[i] + 1]; + m++) { + if (col_P[m] == n) { present = true; - if(n <= col_P[i]) { // make sure we do not add elements twice - sym_col_P[sym_row_P[n] + offset[n]] = col_P[i]; + if (n <= + col_P[i]) { // make sure we do not add elements twice + sym_col_P[sym_row_P[n] + offset[n]] = col_P[i]; sym_col_P[sym_row_P[col_P[i]] + offset[col_P[i]]] = n; - sym_val_P[sym_row_P[n] + offset[n]] = val_P[i] + val_P[m]; - sym_val_P[sym_row_P[col_P[i]] + offset[col_P[i]]] = val_P[i] + val_P[m]; + sym_val_P[sym_row_P[n] + offset[n]] = + val_P[i] + val_P[m]; + sym_val_P[sym_row_P[col_P[i]] + offset[col_P[i]]] = + val_P[i] + val_P[m]; } } } // If (col_P[i], n) is not present, there is no addition involved - if(!present) { - sym_col_P[sym_row_P[n] + offset[n]] = col_P[i]; + if (!present) { + sym_col_P[sym_row_P[n] + offset[n]] = col_P[i]; sym_col_P[sym_row_P[col_P[i]] + offset[col_P[i]]] = n; - sym_val_P[sym_row_P[n] + offset[n]] = val_P[i]; + sym_val_P[sym_row_P[n] + offset[n]] = val_P[i]; sym_val_P[sym_row_P[col_P[i]] + offset[col_P[i]]] = val_P[i]; } // Update offsets - if(!present || (present && n <= col_P[i])) { + if (!present || (present && n <= col_P[i])) { offset[n]++; - if(col_P[i] != n) offset[col_P[i]]++; + if (col_P[i] != n) + offset[col_P[i]]++; } } } // Divide the result by two - for(int i = 0; i < no_elem; i++) sym_val_P[i] /= 2.0; + for (int i = 0; i < no_elem; i++) sym_val_P[i] /= 2.0; // Return symmetrized matrices - free(*_row_P); *_row_P = sym_row_P; - free(*_col_P); *_col_P = sym_col_P; - free(*_val_P); *_val_P = sym_val_P; - - // Free up some memery - free(offset); offset = NULL; - free(row_counts); row_counts = NULL; + *_row_P = move(sym_row_P); + *_col_P = move(sym_col_P); + *_val_P = move(sym_val_P); } // Compute squared Euclidean distance matrix -void TSNE::computeSquaredEuclideanDistance(double* X, int N, int D, double* DD) { +template +void computeSquaredEuclideanDistance(const double* X, int N, double* DD) { const double* XnD = X; - for(int n = 0; n < N; ++n, XnD += D) { + for (int n = 0; n < N; ++n, XnD += D) { const double* XmD = XnD + D; - double* curr_elem = &DD[n*N + n]; + double* curr_elem = &DD[n * N + n]; *curr_elem = 0.0; double* curr_elem_sym = curr_elem + N; - for(int m = n + 1; m < N; ++m, XmD+=D, curr_elem_sym+=N) { + for (int m = n + 1; m < N; ++m, XmD += D, curr_elem_sym += N) { *(++curr_elem) = 0.0; - for(int d = 0; d < D; ++d) { + for (int d = 0; d < D; ++d) { *curr_elem += (XnD[d] - XmD[d]) * (XnD[d] - XmD[d]); } *curr_elem_sym = *curr_elem; @@ -635,156 +672,109 @@ void TSNE::computeSquaredEuclideanDistance(double* X, int N, int D, double* DD) } } +void computeSquaredEuclideanDistance(const double* X, int N, int D, + double* DD) { + const double* XnD = X; + for (int n = 0; n < N; ++n, XnD += D) { + const double* XmD = XnD + D; + double* curr_elem = &DD[n * N + n]; + *curr_elem = 0.0; + double* curr_elem_sym = curr_elem + N; + for (int m = n + 1; m < N; ++m, XmD += D, curr_elem_sym += N) { + *(++curr_elem) = 0.0; + for (int d = 0; d < D; ++d) { + *curr_elem += (XnD[d] - XmD[d]) * (XnD[d] - XmD[d]); + } + *curr_elem_sym = *curr_elem; + } + } +} // Makes data zero-mean -void TSNE::zeroMean(double* X, int N, int D) { - - // Compute data mean - double* mean = (double*) calloc(D, sizeof(double)); - if(mean == NULL) { fprintf(stderr,"Memory allocation failed!\n"); exit(1); } +void zeroMean(double* X, int N, int D) { + // Compute data mean + vector mean(D, 0); int nD = 0; - for(int n = 0; n < N; n++) { - for(int d = 0; d < D; d++) { - mean[d] += X[nD + d]; - } + for (int n = 0; n < N; n++) { + for (int d = 0; d < D; d++) { + mean[d] += X[nD + d]; + } nD += D; - } - for(int d = 0; d < D; d++) { - mean[d] /= (double) N; - } + } + for (int d = 0; d < D; d++) { + mean[d] /= (double)N; + } - // Subtract data mean + // Subtract data mean nD = 0; - for(int n = 0; n < N; n++) { - for(int d = 0; d < D; d++) { - X[nD + d] -= mean[d]; - } + for (int n = 0; n < N; n++) { + for (int d = 0; d < D; d++) { + X[nD + d] -= mean[d]; + } nD += D; - } - free(mean); mean = NULL; -} - - -// Generates a Gaussian random number -double TSNE::randn() { - double x, y, radius; - do { - x = 2 * (rand() / ((double) RAND_MAX + 1)) - 1; - y = 2 * (rand() / ((double) RAND_MAX + 1)) - 1; - radius = (x * x) + (y * y); - } while((radius >= 1.0) || (radius == 0.0)); - radius = sqrt(-2 * log(radius) / radius); - x *= radius; - y *= radius; - return x; -} - -// Function that loads data from a t-SNE file -// Note: this function does a malloc that should be freed elsewhere -bool TSNE::load_data(const char* dat_file, double** data, int* n, int* d, int* no_dims, double* theta, double* perplexity, int* rand_seed, int* max_iter) { - - // Open file, read first 2 integers, allocate memory, and read the data - FILE *h; - if((h = fopen(dat_file, "r+b")) == NULL) { - fprintf(stderr,"Error: could not open data file.\n"); - return false; - } - fread(n, sizeof(int), 1, h); // number of datapoints - fread(d, sizeof(int), 1, h); // original dimensionality - fread(theta, sizeof(double), 1, h); // gradient accuracy - fread(perplexity, sizeof(double), 1, h); // perplexity - fread(no_dims, sizeof(int), 1, h); // output dimensionality - fread(max_iter, sizeof(int),1,h); // maximum number of iterations - *data = (double*) malloc(*d * *n * sizeof(double)); - if(*data == NULL) { fprintf(stderr,"Memory allocation failed!\n"); exit(1); } - fread(*data, sizeof(double), *n * *d, h); // the data - if(!feof(h)) fread(rand_seed, sizeof(int), 1, h); // random seed - fclose(h); - fprintf(stderr,"Read the %i x %i data matrix successfully!\n", *n, *d); - return true; -} - -// Function that saves map to a t-SNE file -void TSNE::save_data(const char* res_file, double* data, int* landmarks, double* costs, int n, int d) { - - // Open file, write first 2 integers and then the data - FILE *h; - if((h = fopen(res_file, "w+b")) == NULL) { - fprintf(stderr,"Error: could not open data file.\n"); - return; - } - fwrite(&n, sizeof(int), 1, h); - fwrite(&d, sizeof(int), 1, h); - fwrite(data, sizeof(double), n * d, h); - fwrite(landmarks, sizeof(int), n, h); - fwrite(costs, sizeof(double), n, h); - fclose(h); - fprintf(stderr,"Wrote the %i x %i data matrix successfully!\n", n, d); + } } -void TSNE::save_csv(const char* csv_file, double* Y, int N, int D) { - std::ofstream csv(csv_file); - +// Makes data zero-mean +template +void zeroMean(double* X, int N) { + // Compute data mean + array mean; + mean.fill(0); + int nD = 0; + for (int n = 0; n < N; n++) { + for (int d = 0; d < D; d++) { + mean[d] += X[nD + d]; + } + nD += D; + } for (int d = 0; d < D; d++) { - csv << "TSNE" << d+1 << ","; + mean[d] /= (double)N; } - csv << "\n"; + // Subtract data mean + nD = 0; for (int n = 0; n < N; n++) { - int row_offset = n * D; for (int d = 0; d < D; d++) { - csv << Y[row_offset + d] << ","; + X[nD + d] -= mean[d]; } - csv << "\n"; + nD += D; } - - csv.close(); } -// Function that runs the Barnes-Hut implementation of t-SNE -int main(int argc, char *argv[]) { +// Generates a Gaussian random number +double randn() { + double x, y, radius; + do { + x = 2 * (rand() / ((double)RAND_MAX + 1)) - 1; + y = 2 * (rand() / ((double)RAND_MAX + 1)) - 1; + radius = (x * x) + (y * y); + } while ((radius >= 1.0) || (radius == 0.0)); + radius = sqrt(-2 * log(radius) / radius); + x *= radius; + y *= radius; + return x; +} - // load input and output - std::string dat_file = "data.dat"; - std::string res_file = "result.dat"; - if (argc > 1) { - dat_file = argv[1]; - res_file = argv[2]; - } +} // namespace - const char *dat_file_c = dat_file.c_str(); - const char *res_file_c = res_file.c_str(); - - // Define some variables - int origN, N, D, no_dims, max_iter, *landmarks; - double perc_landmarks; - double perplexity, theta, *data; - int rand_seed = -1; - TSNE* tsne = new TSNE(); - - // Read the parameters and the dataset - if(tsne->load_data(dat_file_c, &data, &origN, &D, &no_dims, &theta, &perplexity, &rand_seed, &max_iter)) { - - // Make dummy landmarks - N = origN; - int* landmarks = (int*) malloc(N * sizeof(int)); - if(landmarks == NULL) { fprintf(stderr,"Memory allocation failed!\n"); exit(1); } - for(int n = 0; n < N; n++) landmarks[n] = n; - - // Now fire up the SNE implementation - double* Y = (double*) malloc(N * no_dims * sizeof(double)); - double* costs = (double*) calloc(N, sizeof(double)); - if(Y == NULL || costs == NULL) { fprintf(stderr,"Memory allocation failed!\n"); exit(1); } - tsne->run(data, N, D, Y, no_dims, perplexity, theta, rand_seed, false, NULL, false, max_iter); - - // Save the results - tsne->save_data(res_file_c, Y, landmarks, costs, N, no_dims); - - // Clean up the memory - free(data); data = NULL; - free(Y); Y = NULL; - free(costs); costs = NULL; - free(landmarks); landmarks = NULL; +// Perform t-SNE +void run(double* X, int N, int D, double* Y, int no_dims, double perplexity, + double theta, int rand_seed, bool skip_random_init, double* init, + bool use_init, int max_iter, int stop_lying_iter, + int mom_switch_iter) { + switch (no_dims) { + case 2: + run<2>(X, N, D, Y, perplexity, theta, rand_seed, skip_random_init, + init, use_init, max_iter, stop_lying_iter, mom_switch_iter); + return; + case 3: + run<3>(X, N, D, Y, perplexity, theta, rand_seed, skip_random_init, + init, use_init, max_iter, stop_lying_iter, mom_switch_iter); + return; + default: + assert("no_dims must be 2 or 3"); } - delete(tsne); } + +} // namespace TSNE diff --git a/tsne/bh_sne_src/tsne.h b/tsne/bh_sne_src/tsne.h index 93e203e..27aca03 100644 --- a/tsne/bh_sne_src/tsne.h +++ b/tsne/bh_sne_src/tsne.h @@ -30,35 +30,22 @@ * */ - #ifndef TSNE_H #define TSNE_H +namespace TSNE { -static inline double sign(double x) { return (x == .0 ? .0 : (x < .0 ? -1.0 : 1.0)); } - - -class TSNE -{ -public: - void run(double* X, int N, int D, double* Y, int no_dims, double perplexity, double theta, int rand_seed, - bool skip_random_init, double *init, bool use_init, int max_iter=1000, int stop_lying_iter=250, int mom_switch_iter=250 - ); - bool load_data(const char* dat_file, double** data, int* n, int* d, int* no_dims, double* theta, double* perplexity, int* rand_seed, int* max_iter); - void save_data(const char* res_file, double* data, int* landmarks, double* costs, int n, int d); - void symmetrizeMatrix(unsigned int** row_P, unsigned int** col_P, double** val_P, int N); // should be static! - void save_csv(const char* csv_file, double* Y, int N, int D); +void run(double* X, int N, int D, double* Y, int no_dims, double perplexity, + double theta, int rand_seed, bool skip_random_init, double* init, + bool use_init, int max_iter = 1000, int stop_lying_iter = 250, + int mom_switch_iter = 250); +bool load_data(const char* dat_file, double** data, int* n, int* d, + int* no_dims, double* theta, double* perplexity, int* rand_seed, + int* max_iter); +void save_data(const char* res_file, double* data, int* landmarks, + double* costs, int n, int d); +void save_csv(const char* csv_file, double* Y, int N, int D); -private: - 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 computeExactGradient(double* P, double* Y, int N, int D, double* dC); - double evaluateError(double* P, double* Y, int N, int D); - double evaluateError(unsigned int* row_P, unsigned int* col_P, double* val_P, double* Y, int N, int D, double theta); - void zeroMean(double* X, int N, int D); - void computeGaussianPerplexity(double* X, int N, int D, double* P, double perplexity); - void computeGaussianPerplexity(double* X, int N, int D, unsigned int** _row_P, unsigned int** _col_P, double** _val_P, double perplexity, int K); - void computeSquaredEuclideanDistance(double* X, int N, int D, double* DD); - double randn(); -}; +} // namespace TSNE #endif diff --git a/tsne/bh_sne_src/vptree.h b/tsne/bh_sne_src/vptree.h index 9e301d6..7e1f819 100644 --- a/tsne/bh_sne_src/vptree.h +++ b/tsne/bh_sne_src/vptree.h @@ -33,89 +33,60 @@ /* This code was adopted with minor modifications from Steve Hanov's great tutorial at http://stevehanov.ca/blog/index.php?id=130 */ -#include -#include -#include -#include -#include -#include -#include - - #ifndef VPTREE_H #define VPTREE_H -class DataPoint +#include // for reverse, nth_element +#include // for DBL_MAX +#include // for sqrt +#include // for malloc, free, rand, RAND_MAX +#include // for unique_ptr +#include // for priority_queue +#include // for vector + +namespace TSNE { + +class DataPoint final { - int _ind; public: - double* _x; - int _D; - DataPoint() { - _D = 1; - _ind = -1; - _x = NULL; - } - DataPoint(int D, int ind, double* x) { - _D = D; - _ind = ind; - _x = (double*) malloc(_D * sizeof(double)); - for(int d = 0; d < _D; d++) _x[d] = x[d]; - } - DataPoint(const DataPoint& other) { // this makes a deep copy -- should not free anything - if(this != &other) { - _D = other.dimensionality(); - _ind = other.index(); - _x = (double*) malloc(_D * sizeof(double)); - for(int d = 0; d < _D; d++) _x[d] = other.x(d); - } - } - ~DataPoint() { if(_x != NULL) free(_x); } - DataPoint& operator= (const DataPoint& other) { // asignment should free old object - if(this != &other) { - if(_x != NULL) free(_x); - _D = other.dimensionality(); - _ind = other.index(); - _x = (double*) malloc(_D * sizeof(double)); - for(int d = 0; d < _D; d++) _x[d] = other.x(d); - } - return *this; - } - int index() const { return _ind; } - int dimensionality() const { return _D; } + const double* _x; + DataPoint() = default; + DataPoint(const double* x) : _x(x) {} + DataPoint(const DataPoint& other) = default; double x(int d) const { return _x[d]; } }; -double euclidean_distance(const DataPoint &t1, const DataPoint &t2) { - double dd = .0; - double* x1 = t1._x; - double* x2 = t2._x; - double diff; - for(int d = 0; d < t1._D; d++) { - diff = (x1[d] - x2[d]); - dd += diff * diff; +class euclidean_distance final { + const int _D; + public: + euclidean_distance(int D) : _D(D) {} + euclidean_distance(euclidean_distance&&) = default; + euclidean_distance(const euclidean_distance&) = default; + double operator()(const DataPoint &t1, const DataPoint &t2) const { + double dd = .0; + const double * x1 = t1._x; + const double * x2 = t2._x; + double diff; + for(int d = 0; d < _D; d++) { + diff = (x1[d] - x2[d]); + dd += diff * diff; + } + return sqrt(dd); } - return sqrt(dd); -} +}; -template -class VpTree +template +class VpTree final { + VpTree(const VpTree&) = delete; public: - // Default constructor - VpTree() : _root(0) {} + explicit VpTree(Distance&& distance) : distance(distance) {} - // Destructor - ~VpTree() { - delete _root; - } - // Function to create a new VpTree from data void create(const std::vector& items) { - delete _root; _items = items; _root = buildFromPoints(0, items.size()); } @@ -131,7 +102,7 @@ class VpTree _tau = DBL_MAX; // Perform the search - search(_root, target, k, heap); + search(_root.get(), target, k, heap); // Gather final results results->clear(); distances->clear(); @@ -147,6 +118,7 @@ class VpTree } private: + const Distance distance; std::vector _items; double _tau; @@ -155,17 +127,12 @@ class VpTree { int index; // index of point in node double threshold; // radius(?) - Node* left; // points closer by than threshold - Node* right; // points farther away than threshold + std::unique_ptr left; // points closer by than threshold + std::unique_ptr right; // points farther away than threshold - Node() : - index(0), threshold(0.), left(0), right(0) {} - - ~Node() { // destructor - delete left; - delete right; - } - }* _root; + Node() : index(0), threshold(0.) {} + }; + std::unique_ptr _root; // An item on the intermediate result queue @@ -183,21 +150,23 @@ class VpTree struct DistanceComparator { const T& item; - DistanceComparator(const T& item) : item(item) {} - bool operator()(const T& a, const T& b) { + const Distance distance; + DistanceComparator(const T& item, + const Distance distance) : item(item), distance(distance) {} + bool operator()(const T& a, const T& b) const { return distance(item, a) < distance(item, b); } }; // Function that (recursively) fills the tree - Node* buildFromPoints( int lower, int upper ) + std::unique_ptr buildFromPoints( int lower, int upper ) { if (upper == lower) { // indicates that we're done here! - return NULL; + return std::unique_ptr(); } // Lower index is center of current node - Node* node = new Node(); + std::unique_ptr node(new Node()); node->index = lower; if (upper - lower > 1) { // if we did not arrive at leaf yet @@ -211,7 +180,7 @@ class VpTree std::nth_element(_items.begin() + lower + 1, _items.begin() + median, _items.begin() + upper, - DistanceComparator(_items[lower])); + DistanceComparator(_items[lower], distance)); // Threshold of the new node will be the distance to the median node->threshold = distance(_items[lower], _items[median]); @@ -229,7 +198,7 @@ class VpTree // Helper function that searches the tree void search(Node* node, const T& target, int k, std::priority_queue& heap) { - if(node == NULL) return; // indicates that we're done here + if(node == nullptr) return; // indicates that we're done here // Compute distance between target and current node double dist = distance(_items[node->index], target); @@ -242,31 +211,33 @@ class VpTree } // Return if we arrived at a leaf - if(node->left == NULL && node->right == NULL) { + if(node->left == nullptr && node->right == nullptr) { return; } // If the target lies within the radius of ball if(dist < node->threshold) { if(dist - _tau <= node->threshold) { // if there can still be neighbors inside the ball, recursively search left child first - search(node->left, target, k, heap); + search(node->left.get(), target, k, heap); } if(dist + _tau >= node->threshold) { // if there can still be neighbors outside the ball, recursively search right child - search(node->right, target, k, heap); + search(node->right.get(), target, k, heap); } // If the target lies outsize the radius of the ball } else { if(dist + _tau >= node->threshold) { // if there can still be neighbors outside the ball, recursively search right child first - search(node->right, target, k, heap); + search(node->right.get(), target, k, heap); } if (dist - _tau <= node->threshold) { // if there can still be neighbors inside the ball, recursively search left child - search(node->left, target, k, heap); + search(node->left.get(), target, k, heap); } } } }; + +} // namespace TSNE #endif