From 5347229be3a9431b8a8d1a3ee75f37314cee454b Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Thu, 31 May 2018 23:38:37 -0700 Subject: [PATCH 01/22] Don't override CXX, CFLAGS from environment. Also, enable link-time optimize. --- tsne/bh_sne_src/Makefile | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/tsne/bh_sne_src/Makefile b/tsne/bh_sne_src/Makefile index 8c8a1c8..20b37f2 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 + +CFLAGS += -ffunction-sections +LDFLAGS += -Wl,--gc-sections all: bh_tsne bh_tsne_3d bh_tsne: tsne.o sptree.o - $(CXX) $(CFLAGS) tsne.o sptree.o -o bh_tsne + $(CXX) $(CFLAGS) $(LDFLAGS) $^ -o $@ bh_tsne_3d: tsne_3d.o sptree.o - $(CXX) $(CFLAGS) tsne_3d.o sptree.o -o bh_tsne_3d + $(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 + $(CXX) $(CFLAGS) -DTSNE3D -c $< -o $@ clean: rm -Rf *.o bh_tsne bh_tsne_3d From c0de07f4a70c3a3f8fb7b9825b0d32aac7fc68e2 Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Thu, 31 May 2018 23:45:37 -0700 Subject: [PATCH 02/22] Optimize and clean up SPTree. Inspired by https://www.microsoft.com/en-us/research/blog/optimizing-barnes-hut-t-sne/ but backported to our fork. First of all, remove all manual memory management and use STL containers as appropriate. Remove the buffer used for caching diffs, and instead compute on the fly - it's a cheap computation, and repeating it reduces memory cache misses which matters a lot more than an extra subtraction or two. It also saves doubles worth of memory per node. Share the widths structure between nodes at the same level. This saves ( - 2) * 2^ pointers, and costs doubles. So at dimension 2 this actually costs more than it saves, but still helps with memory locality. At higher dimension it saves a lot. Minor optimization of the forces inequality, which should help both memory locality and speed. Ran include-what-you-use to get the headers right. Removed using namespace std to avoid problems with future C++ standards upgrades. --- tsne/bh_sne_src/sptree.cpp | 243 +++++++++++++++++++------------------ tsne/bh_sne_src/sptree.h | 85 ++++++++----- 2 files changed, 178 insertions(+), 150 deletions(-) diff --git a/tsne/bh_sne_src/sptree.cpp b/tsne/bh_sne_src/sptree.cpp index 2d187c6..40b03e2 100644 --- a/tsne/bh_sne_src/sptree.cpp +++ b/tsne/bh_sne_src/sptree.cpp @@ -30,78 +30,69 @@ * */ -#include -#include -#include -#include -#include #include "sptree.h" +#include // for move, max +#include // for DBL_MAX +#include // for fprintf, stderr, size_t +#include // for unique_ptr +using std::max; +using std::max_element; +using std::move; +using std::unique_ptr; // Constructs cell template -Cell::Cell() { -} - -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]); -} - -// Destructs cell -template -Cell::~Cell() { -} +Cell::Cell(typename Cell::point_t&& inp_corner, + const typename Cell::point_t* width) : + width(width), corner(inp_corner) {} template -double Cell::getCorner(unsigned int d) { +double Cell::getCorner(unsigned int d) const { return corner[d]; } template -double Cell::getWidth(unsigned int d) { - return width[d]; +double Cell::getWidth(unsigned int d) const { + return (*width)[d]; } template -void Cell::setCorner(unsigned int d, double val) { - corner[d] = val; +void Cell::setCorner(typename Cell::point_t&& val) { + corner = val; } template -void Cell::setWidth(unsigned int d, double val) { - width[d] = val; +void Cell::setWidth(const typename Cell::point_t* val) { + width = val; } // Checks whether a point lies in a cell template -bool Cell::containsPoint(double point[]) +bool Cell::containsPoint(const double* point) const { 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; + if(corner[d] - (*width)[d] > point[d]) return false; + if(corner[d] + (*width)[d] < point[d]) return false; } return true; } - -// Default constructor for SPTree -- build tree, too! template -SPTree::SPTree(double* inp_data, unsigned int N) -{ +double Cell::maxWidth() const { + return *max_element(width->begin(), width->end()); +} +// Top-node constructor for SPTree -- build tree, too! +template +SPTree::SPTree(const 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; - } + point_t mean_Y, min_Y, max_Y; + mean_Y.fill(0.0); + min_Y.fill(DBL_MAX); + max_Y.fill(-DBL_MAX); for(unsigned int n = 0; n < N; n++) { for(unsigned int d = 0; d < NDims; d++) { @@ -115,55 +106,63 @@ SPTree::SPTree(double* inp_data, unsigned int N) 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); + top_widths = unique_ptr(new point_t()); + auto& widths = *top_widths; + for(int d = 0; d < NDims; d++) { + widths[d] = max(max_Y[d] - mean_Y[d], mean_Y[d] - min_Y[d]) + 1e-5; + } + init(nullptr, inp_data, move(mean_Y), top_widths.get()); fill(N); - - // Clean up memory - free(mean_Y); - free(max_Y); - free(min_Y); - free(width); } // 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) +SPTree::SPTree(const double* inp_data, + unsigned int N, + typename SPTree::point_t&& inp_corner, + const typename SPTree::point_t* inp_width) { - init(NULL, inp_data, inp_corner, inp_width); + init(nullptr, inp_data, move(inp_corner), inp_width); fill(N); } // Constructor for SPTree with particular size (do not fill the tree) template -SPTree::SPTree(double* inp_data, double* inp_corner, double* inp_width) +SPTree::SPTree(const double* inp_data, + typename SPTree::point_t&& inp_corner, + const typename SPTree::point_t* inp_width) { - init(NULL, inp_data, inp_corner, inp_width); + init(nullptr, inp_data, move(inp_corner), inp_width); } // 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); +SPTree::SPTree(SPTree* inp_parent, const double* inp_data, + typename SPTree::point_t&& inp_corner, + const typename SPTree::point_t* inp_width) { + init(inp_parent, inp_data, move(inp_corner), inp_width); } // 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) +SPTree::SPTree(SPTree* inp_parent, const double* inp_data, unsigned int N, + typename SPTree::point_t&& inp_corner, + const typename SPTree::point_t* inp_width) { - init(inp_parent, inp_data, inp_corner, inp_width); + init(inp_parent, inp_data, move(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) +void SPTree::init(SPTree* inp_parent, const double* inp_data, + typename SPTree::point_t&& inp_corner, + const typename SPTree::point_t* inp_width) { parent = inp_parent; data = inp_data; @@ -171,27 +170,16 @@ void SPTree::init(SPTree* inp_parent, double* inp_data, double* inp_corne 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]); + boundary.setCorner(move(inp_corner)); + boundary.setWidth(inp_width); - 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]; - } + center_of_mass.fill(0.0); } // Update the data underlying this tree template -void SPTree::setData(double* inp_data) +void SPTree::setData(const double* inp_data) { data = inp_data; } @@ -210,7 +198,7 @@ template bool SPTree::insert(unsigned int new_index) { // Ignore objects which do not belong in this quad tree - double* point = data + new_index * NDims; + const double* point = data + new_index * NDims; if(!boundary.containsPoint(point)) return false; @@ -245,8 +233,8 @@ bool SPTree::insert(unsigned int new_index) 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; + for(auto& child : children) { + if(child->insert(new_index)) return true; } // Otherwise, the point cannot be inserted (this should never happen) @@ -259,26 +247,32 @@ template void SPTree::subdivide() { // Create new children - double new_corner[NDims]; - double new_width[NDims]; + child_widths = unique_ptr(new point_t()); + point_t& new_width = *child_widths; + for(unsigned int d = 0; d < NDims; d++) { + new_width[d] = .5 * boundary.getWidth(d); + } + for(unsigned int i = 0; i < no_children; i++) { unsigned int div = 1; + point_t new_corner; 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); + if((i / div) % 2 == 1) + new_corner[d] = boundary.getCorner(d) - new_width[d]; + else + new_corner[d] = boundary.getCorner(d) + new_width[d]; div *= 2; } - children[i] = new SPTree(this, data, new_corner, new_width); + children[i] = unique_ptr(new SPTree(this, data, move(new_corner), child_widths.get())); } // 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]); + for (auto& child : children) { + if (child->insert(index[i])) { + break; + } } - index[i] = -1; } // Empty parent node @@ -297,15 +291,15 @@ void SPTree::fill(unsigned int N) // Checks whether the specified tree is correct template -bool SPTree::isCorrect() +bool SPTree::isCorrect() const { for(unsigned int n = 0; n < size; n++) { - double* point = data + index[n] * NDims; + const 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(); + for(const auto& child : children) correct = correct && child->isCorrect(); return correct; } else return true; @@ -315,7 +309,7 @@ bool SPTree::isCorrect() // Build a list of all indices in SPTree template -void SPTree::getAllIndices(unsigned int* indices) +void SPTree::getAllIndices(unsigned int* indices) const { getAllIndices(indices, 0); } @@ -323,7 +317,7 @@ void SPTree::getAllIndices(unsigned int* indices) // Build a list of all indices in SPTree template -unsigned int SPTree::getAllIndices(unsigned int* indices, unsigned int loc) +unsigned int SPTree::getAllIndices(unsigned int* indices, unsigned int loc) const { // Gather indices in current quadrant @@ -332,23 +326,29 @@ unsigned int SPTree::getAllIndices(unsigned int* indices, unsigned int lo // Gather indices in children if(!is_leaf) { - for(int i = 0; i < no_children; i++) loc = children[i]->getAllIndices(indices, loc); + for(const auto& child : children) loc = child->getAllIndices(indices, loc); } return loc; } template -unsigned int SPTree::getDepth() { +unsigned int SPTree::getDepth() const { 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; + 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) { + double max_width = boundary.maxWidth(); + computeNonEdgeForces(point_index, theta, neg_f, sum_Q, max_width*max_width); +} // Compute non-edge forces using Barnes-Hut algorithm template -void SPTree::computeNonEdgeForces(unsigned int point_index, double theta, double neg_f[], double* sum_Q) +void SPTree::computeNonEdgeForces(unsigned int point_index, double theta, double neg_f[], double* sum_Q, double max_width_squared) { // Make sure that we spend no time on empty nodes or self-interactions @@ -356,33 +356,40 @@ void SPTree::computeNonEdgeForces(unsigned int point_index, double theta, // Compute distance between point and center-of-mass double sqdist = .0; - unsigned int ind = point_index * NDims; + size_t 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(const auto& cm : center_of_mass) { + double diff = data[ind] - cm; + ind += NDims; + sqdist += diff * diff; } // 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) { + // max_width / sqrt(sqdist) < theta + if(is_leaf || max_width_squared < theta * theta * sqdist) { // 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]; + ind = point_index * NDims; + for(size_t d = 0; d < NDims; ++d) { + // recompute here rather than storing from before because memory + // locality matters more than an extra couple of additions and a + // subtraction. + double diff = data[ind] - center_of_mass[d]; + ind += NDims; + neg_f[d] += mult * diff; + } } 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); + max_width_squared /= 4.0; + for(unsigned int i = 0; i < no_children; i++) + children[i]->computeNonEdgeForces( + point_index, theta, neg_f, sum_Q, max_width_squared); } } @@ -394,24 +401,26 @@ void SPTree::computeEdgeForces(unsigned int* row_P, unsigned int* col_P, // 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++) { // Compute pairwise distance and Q-value sqdist = 1.0; - ind2 = col_P[i] * NDims; + unsigned int ind2 = 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]; + double diff = data[ind1 + d] - data[ind2 + d]; + sqdist += diff * diff; } 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++) { + double diff = data[ind1 + d] - data[ind2 + d]; + pos_f[ind1 + d] += sqdist * diff; + } } ind1 += NDims; } @@ -430,7 +439,7 @@ void SPTree::print() if(is_leaf) { fprintf(stderr,"Leaf node; data = ["); for(int i = 0; i < size; i++) { - double* point = data + index[i] * NDims; + 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 < size - 1) fprintf(stderr,"\n"); @@ -439,7 +448,7 @@ void SPTree::print() } else { fprintf(stderr,"Intersection node with center-of-mass = ["); - for(int d = 0; d < NDims; d++) fprintf(stderr,"%f, ", center_of_mass[d]); + for(const auto& cm : center_of_mass) fprintf(stderr,"%f, ", cm); fprintf(stderr,"]; children are:\n"); for(int i = 0; i < no_children; i++) children[i]->print(); } diff --git a/tsne/bh_sne_src/sptree.h b/tsne/bh_sne_src/sptree.h index be42686..9380774 100644 --- a/tsne/bh_sne_src/sptree.h +++ b/tsne/bh_sne_src/sptree.h @@ -34,25 +34,30 @@ #ifndef SPTREE_H #define SPTREE_H -using namespace std; - +#include +#include template class Cell { - double corner[NDims]; - double width[NDims]; + public: + typedef std::array point_t; + Cell() = default; + Cell(point_t&& inp_corner, const point_t* width); -public: - Cell(); - Cell(double* inp_corner, double* inp_width); - ~Cell(); - - 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[]); + double getCorner(unsigned int d) const; + double getWidth(unsigned int d) const; + void setCorner(point_t&& inp_corner); + void setWidth(const point_t* width); + bool containsPoint(const double* point) const; + double maxWidth() const; + + private: + point_t corner; + const point_t* width; + + // disallow copy + Cell(const Cell&) = delete; }; template @@ -61,13 +66,12 @@ class SPTree public: enum { no_children = 2 * SPTree::no_children }; + typedef typename Cell::point_t point_t; + 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; @@ -75,41 +79,56 @@ class SPTree unsigned int size; unsigned int cum_size; + // The width for each cell is the same at each level. The parent + // node owns the widths array for its children. The top node must + // also own the array for itself. + std::unique_ptr top_widths; + std::unique_ptr child_widths; + // 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]; + const double* data; + point_t center_of_mass; + std::array index; // Children - SPTree* children[no_children]; + std::array>, no_children> children; + + SPTree(const double* inp_data, point_t&& inp_corner, const point_t* inp_width); + SPTree(const double* inp_data, unsigned int N, point_t&& inp_corner, const point_t* inp_width); + SPTree(SPTree* inp_parent, const double* inp_data, unsigned int N, point_t&& inp_corner, const point_t* inp_width); + SPTree(SPTree* inp_parent, const double* inp_data, point_t&& inp_corner, const point_t* inp_width); + + // Disallow copy + SPTree(const SPTree&) = delete; 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(const double* inp_data, unsigned int N); + + void setData(const double* inp_data); SPTree* getParent(); void construct(Cell boundary); bool insert(unsigned int new_index); void subdivide(); - bool isCorrect(); + bool isCorrect() const; void rebuildTree(); - void getAllIndices(unsigned int* indices); - unsigned int getDepth(); + void getAllIndices(unsigned int* indices) const; + unsigned int getDepth() const; 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 init(SPTree* inp_parent, const double* inp_data, point_t&& inp_corner, const point_t* inp_width); + void computeNonEdgeForces(unsigned int point_index, + double theta, + double neg_f[], + double* sum_Q, + double max_width_squared); void fill(unsigned int N); - unsigned int getAllIndices(unsigned int* indices, unsigned int loc); + unsigned int getAllIndices(unsigned int* indices, unsigned int loc) const; bool isChild(unsigned int test_index, unsigned int start, unsigned int end); }; From 94be39c93d3c5e88e6e6a9d045ec69a33affa43a Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Fri, 1 Jun 2018 02:44:19 -0700 Subject: [PATCH 03/22] Stop carrying around top-node stuff in child nodes. Also stop carrying around parent pointers. We don't need them. Stuff that can be easily flowed down from the top at query time now is, rather than being duplicated at every level. All told, by special-casing the top node we can remove 4 pointers, an int, and a bool from every node. --- tsne/bh_sne_src/sptree.cpp | 239 ++++++++++++++++--------------------- tsne/bh_sne_src/sptree.h | 116 ++++++++++-------- 2 files changed, 166 insertions(+), 189 deletions(-) diff --git a/tsne/bh_sne_src/sptree.cpp b/tsne/bh_sne_src/sptree.cpp index 40b03e2..6fe8dd6 100644 --- a/tsne/bh_sne_src/sptree.cpp +++ b/tsne/bh_sne_src/sptree.cpp @@ -32,56 +32,45 @@ #include "sptree.h" -#include // for move, max +#include // for max, max_element #include // for DBL_MAX #include // for fprintf, stderr, size_t #include // for unique_ptr +#include // for move using std::max; using std::max_element; using std::move; using std::unique_ptr; +using std::vector; -// Constructs cell template -Cell::Cell(typename Cell::point_t&& inp_corner, - const typename Cell::point_t* width) : - width(width), corner(inp_corner) {} +Cell::Cell(typename Cell::point_t&& p) : corner(move(p)) {} template double Cell::getCorner(unsigned int d) const { return corner[d]; } -template -double Cell::getWidth(unsigned int d) const { - return (*width)[d]; -} - template void Cell::setCorner(typename Cell::point_t&& val) { corner = val; } -template -void Cell::setWidth(const typename Cell::point_t* val) { - width = val; -} - // Checks whether a point lies in a cell template -bool Cell::containsPoint(const double* point) const +bool Cell::containsPoint(const double* point, const typename Cell::point_t& width) const { 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; + if(corner[d] - width[d] > point[d]) return false; + if(corner[d] + width[d] < point[d]) return false; } return true; } template -double Cell::maxWidth() const { - return *max_element(width->begin(), width->end()); +double SPTree::maxWidth() const { + return *max_element(widths[0].begin(), widths[0].end()); } // Top-node constructor for SPTree -- build tree, too! @@ -103,103 +92,53 @@ SPTree::SPTree(const double* inp_data, unsigned int N) { nD += NDims; } - for(int d = 0; d < NDims; d++) mean_Y[d] /= (double) N; + double dbl_N = static_cast(N); + for(int d = 0; d < NDims; d++) mean_Y[d] /= dbl_N; // Construct SPTree - top_widths = unique_ptr(new point_t()); - auto& widths = *top_widths; + point_t width; for(int d = 0; d < NDims; d++) { - widths[d] = max(max_Y[d] - mean_Y[d], mean_Y[d] - min_Y[d]) + 1e-5; + width[d] = max(max_Y[d] - mean_Y[d], mean_Y[d] - min_Y[d]) + 1e-5; } - init(nullptr, inp_data, move(mean_Y), top_widths.get()); - fill(N); -} - - -// Constructor for SPTree with particular size and parent -- build the tree, too! -template -SPTree::SPTree(const double* inp_data, - unsigned int N, - typename SPTree::point_t&& inp_corner, - const typename SPTree::point_t* inp_width) -{ - init(nullptr, inp_data, move(inp_corner), inp_width); - fill(N); -} - - -// Constructor for SPTree with particular size (do not fill the tree) -template -SPTree::SPTree(const double* inp_data, - typename SPTree::point_t&& inp_corner, - const typename SPTree::point_t* inp_width) -{ - init(nullptr, inp_data, move(inp_corner), inp_width); -} - - -// Constructor for SPTree with particular size and parent (do not fill tree) -template -SPTree::SPTree(SPTree* inp_parent, const double* inp_data, - typename SPTree::point_t&& inp_corner, - const typename SPTree::point_t* inp_width) { - init(inp_parent, inp_data, move(inp_corner), inp_width); -} - - -// Constructor for SPTree with particular size and parent -- build the tree, too! -template -SPTree::SPTree(SPTree* inp_parent, const double* inp_data, unsigned int N, - typename SPTree::point_t&& inp_corner, - const typename SPTree::point_t* inp_width) -{ - init(inp_parent, inp_data, move(inp_corner), inp_width); + widths.emplace_back(move(width)); + node.setCorner(move(mean_Y)); fill(N); } - -// Main initialization function +// Constructor for SPTreeNode. template -void SPTree::init(SPTree* inp_parent, const double* inp_data, - typename SPTree::point_t&& inp_corner, - const typename SPTree::point_t* inp_width) -{ - parent = inp_parent; - data = inp_data; - is_leaf = true; +SPTreeNode::SPTreeNode(typename SPTreeNode::point_t&& inp_corner) : + boundary(move(inp_corner)) { size = 0; cum_size = 0; - boundary.setCorner(move(inp_corner)); - boundary.setWidth(inp_width); - center_of_mass.fill(0.0); } - -// Update the data underlying this tree +// Update the corner position. template -void SPTree::setData(const double* inp_data) -{ - data = inp_data; +void SPTreeNode::setCorner(typename SPTreeNode::point_t&& inp_corner) { + boundary.setCorner(move(inp_corner)); } +template +bool SPTreeNode::is_leaf() const { + return NDims == 0 || children[0].get() == nullptr; +} -// Get the parent of the current tree +// Insert a point into the SPTree template -SPTree* SPTree::getParent() -{ - return parent; +bool SPTree::insert(unsigned int new_index) { + return node.insert(new_index, data, &widths, 1); } - -// Insert a point into the SPTree template -bool SPTree::insert(unsigned int new_index) +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 const double* point = data + new_index * NDims; - if(!boundary.containsPoint(point)) + const auto& width = (*widths)[depth]; + if(!boundary.containsPoint(point, width)) return false; // Online update of cumulative size and center-of-mass @@ -212,29 +151,31 @@ bool SPTree::insert(unsigned int new_index) } // 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 && is_leaf()) { 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; } } - any_duplicate = any_duplicate | duplicate; + if (duplicate) { + return true; + } } - if(any_duplicate) return true; // Otherwise, we need to subdivide the current cell - if(is_leaf) subdivide(); + if(is_leaf()) subdivide(data, widths, depth); // Find out where the point can be inserted for(auto& child : children) { - if(child->insert(new_index)) return true; + if(child->insert(new_index, data, widths, depth+1)) { + return true; + } } // Otherwise, the point cannot be inserted (this should never happen) @@ -244,15 +185,21 @@ bool SPTree::insert(unsigned int new_index) // Create four children which fully divide this cell into four quads of equal area template -void SPTree::subdivide() { +void SPTreeNode::subdivide(const double* data, vector* widths, typename vector::size_type depth) { - // Create new children - child_widths = unique_ptr(new point_t()); - point_t& new_width = *child_widths; - for(unsigned int d = 0; d < NDims; d++) { - new_width[d] = .5 * boundary.getWidth(d); + // If nessessary, add to the width. + if (depth+1 == widths->size()) { + // extend the list. + point_t child_widths; + const point_t& width = (*widths)[depth]; + for(unsigned int d = 0; d < NDims; d++) { + child_widths[d] = .5 * width[d]; + } + widths->emplace_back(move(child_widths)); } + const point_t& new_width = (*widths)[depth+1]; + // Create new children for(unsigned int i = 0; i < no_children; i++) { unsigned int div = 1; point_t new_corner; @@ -263,13 +210,13 @@ void SPTree::subdivide() { new_corner[d] = boundary.getCorner(d) + new_width[d]; div *= 2; } - children[i] = unique_ptr(new SPTree(this, data, move(new_corner), child_widths.get())); + children[i].reset(new SPTreeNode(move(new_corner))); } // Move existing points to correct children for(unsigned int i = 0; i < size; i++) { for (auto& child : children) { - if (child->insert(index[i])) { + if (child->insert(index[i], data, widths, depth+1)) { break; } } @@ -277,7 +224,6 @@ void SPTree::subdivide() { // Empty parent node size = 0; - is_leaf = false; } @@ -293,16 +239,26 @@ void SPTree::fill(unsigned int N) template bool SPTree::isCorrect() const { + return node.isCorrect(data, widths.begin()); +} + + +template +bool SPTreeNode::isCorrect(const double* data, + typename vector::const_iterator width) const { for(unsigned int n = 0; n < size; n++) { const double* point = data + index[n] * NDims; - if(!boundary.containsPoint(point)) return false; + if(!boundary.containsPoint(point, *width)) return false; } - if(!is_leaf) { - bool correct = true; - for(const auto& child : children) correct = correct && child->isCorrect(); - return correct; + if(!is_leaf()) { + ++width; + for(const auto& child : children) { + if (!child->isCorrect(data, width)) { + return false; + } + } } - else return true; + return true; } @@ -311,13 +267,13 @@ bool SPTree::isCorrect() const template void SPTree::getAllIndices(unsigned int* indices) const { - getAllIndices(indices, 0); + node.getAllIndices(indices, 0); } // Build a list of all indices in SPTree template -unsigned int SPTree::getAllIndices(unsigned int* indices, unsigned int loc) const +unsigned int SPTreeNode::getAllIndices(unsigned int* indices, unsigned int loc) const { // Gather indices in current quadrant @@ -325,15 +281,21 @@ unsigned int SPTree::getAllIndices(unsigned int* indices, unsigned int lo loc += size; // Gather indices in children - if(!is_leaf) { - for(const auto& child : children) loc = child->getAllIndices(indices, loc); + if(!is_leaf()) { + for(const auto& child : children) + loc = child->getAllIndices(indices, loc); } return loc; } template unsigned int SPTree::getDepth() const { - if(is_leaf) return 1; + return node.getDepth(); +} + +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; @@ -341,45 +303,43 @@ unsigned int SPTree::getDepth() const { // Compute non-edge forces using Barnes-Hut algorithm template -void SPTree::computeNonEdgeForces(unsigned int point_index, double theta, double neg_f[], double* sum_Q) { - double max_width = boundary.maxWidth(); - computeNonEdgeForces(point_index, theta, neg_f, sum_Q, max_width*max_width); +void SPTree::computeNonEdgeForces(unsigned int point_index, double theta, double neg_f[], double* sum_Q) const { + double max_width = maxWidth(); + node.computeNonEdgeForces(point_index, data + point_index * NDims, theta, neg_f, sum_Q, max_width*max_width); } // Compute non-edge forces using Barnes-Hut algorithm template -void SPTree::computeNonEdgeForces(unsigned int point_index, double theta, double neg_f[], double* sum_Q, double max_width_squared) +void SPTreeNode::computeNonEdgeForces(unsigned int point_index, + const double* data_point, double theta, 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 || (is_leaf() && size == 1 && index[0] == point_index)) return; // Compute distance between point and center-of-mass double sqdist = .0; - size_t ind = point_index * NDims; - for(const auto& cm : center_of_mass) { - double diff = data[ind] - cm; - ind += NDims; + for(int i = 0; i < NDims; ++i) { + double diff = data_point[i] - center_of_mass[i]; sqdist += diff * diff; } // Check whether we can use this node as a "summary" // max_width / sqrt(sqdist) < theta - if(is_leaf || max_width_squared < theta * theta * sqdist) { + if(is_leaf() || max_width_squared < theta * theta * sqdist) { // 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; - ind = point_index * NDims; for(size_t d = 0; d < NDims; ++d) { // recompute here rather than storing from before because memory // locality matters more than an extra couple of additions and a // subtraction. - double diff = data[ind] - center_of_mass[d]; - ind += NDims; + double diff = data_point[d] - center_of_mass[d]; neg_f[d] += mult * diff; } } @@ -388,15 +348,15 @@ void SPTree::computeNonEdgeForces(unsigned int point_index, double theta, // Recursively apply Barnes-Hut to children max_width_squared /= 4.0; for(unsigned int i = 0; i < no_children; i++) - children[i]->computeNonEdgeForces( - point_index, theta, neg_f, sum_Q, max_width_squared); + children[i]->computeNonEdgeForces(point_index, data_point, + theta, neg_f, sum_Q, max_width_squared); } } // Computes edge forces template -void SPTree::computeEdgeForces(unsigned int* row_P, unsigned int* col_P, double* val_P, int N, double* pos_f) +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 @@ -426,17 +386,22 @@ void SPTree::computeEdgeForces(unsigned int* row_P, unsigned int* col_P, } } +// Print out tree +template +void SPTree::print() const { + node.print(data); +} // Print out tree template -void SPTree::print() +void SPTreeNode::print(const double* data) const { if(cum_size == 0) { fprintf(stderr,"Empty node\n"); return; } - if(is_leaf) { + if(is_leaf()) { fprintf(stderr,"Leaf node; data = ["); for(int i = 0; i < size; i++) { const double* point = data + index[i] * NDims; @@ -450,7 +415,7 @@ void SPTree::print() 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(int i = 0; i < no_children; i++) children[i]->print(); + for(int i = 0; i < no_children; i++) children[i]->print(data); } } diff --git a/tsne/bh_sne_src/sptree.h b/tsne/bh_sne_src/sptree.h index 9380774..3e7a459 100644 --- a/tsne/bh_sne_src/sptree.h +++ b/tsne/bh_sne_src/sptree.h @@ -36,106 +36,118 @@ #include #include +#include template -class Cell { +class Cell final { public: typedef std::array point_t; Cell() = default; - Cell(point_t&& inp_corner, const point_t* width); + Cell(point_t&&); double getCorner(unsigned int d) const; - double getWidth(unsigned int d) const; void setCorner(point_t&& inp_corner); - void setWidth(const point_t* width); - bool containsPoint(const double* point) const; - double maxWidth() const; + bool containsPoint(const double* point, const point_t& width) const; private: point_t corner; - const point_t* width; // disallow copy Cell(const Cell&) = delete; }; -template -class SPTree -{ +template +class SPTreeNode final { public: - enum { no_children = 2 * SPTree::no_children }; - typedef typename Cell::point_t point_t; + typedef typename Cell::point_t point_t; + enum { no_children = 2 * SPTreeNode::no_children }; private: + // Fixed constants - static const unsigned int QT_NODE_CAPACITY = 1; + 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 + Cell boundary; // Properties of this node in the tree - SPTree* parent; - unsigned int dimension; - bool is_leaf; unsigned int size; unsigned int cum_size; - // The width for each cell is the same at each level. The parent - // node owns the widths array for its children. The top node must - // also own the array for itself. - std::unique_ptr top_widths; - std::unique_ptr child_widths; - - // 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 - const double* data; point_t center_of_mass; std::array index; // Children - std::array>, no_children> children; - - SPTree(const double* inp_data, point_t&& inp_corner, const point_t* inp_width); - SPTree(const double* inp_data, unsigned int N, point_t&& inp_corner, const point_t* inp_width); - SPTree(SPTree* inp_parent, const double* inp_data, unsigned int N, point_t&& inp_corner, const point_t* inp_width); - SPTree(SPTree* inp_parent, const double* inp_data, point_t&& inp_corner, const point_t* inp_width); + std::array>, no_children> children; // Disallow copy - SPTree(const SPTree&) = delete; + SPTreeNode(const SPTreeNode&) = delete; + + void subdivide(const double* data, std::vector* widths, typename std::vector::size_type depth); + + bool is_leaf() const; public: - SPTree(const double* inp_data, unsigned int N); + SPTreeNode() = default; + SPTreeNode(point_t&& corner); - void setData(const double* inp_data); - SPTree* getParent(); - void construct(Cell boundary); - bool insert(unsigned int new_index); - void subdivide(); - bool isCorrect() const; - void rebuildTree(); - void getAllIndices(unsigned int* indices) const; - unsigned int getDepth() const; - 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(); + void setCorner(point_t&& corner); + + bool insert(unsigned int new_index, const double* data, std::vector* widths, typename std::vector::size_type depth); + bool isCorrect(const double* data, typename std::vector::const_iterator width) const; -private: - void init(SPTree* inp_parent, const double* inp_data, point_t&& inp_corner, const point_t* inp_width); void computeNonEdgeForces(unsigned int point_index, + const double* data_point, double theta, double neg_f[], double* sum_Q, - double max_width_squared); - void fill(unsigned int N); + double max_width_squared) const; unsigned int getAllIndices(unsigned int* indices, unsigned int loc) const; - bool isChild(unsigned int test_index, unsigned int start, unsigned int end); + unsigned int getDepth() const; + void print(const double* data) const; }; template <> -struct SPTree<0> +struct SPTreeNode<0> { enum { no_children = 1 }; }; +template +class SPTree +{ +public: + typedef typename SPTreeNode::point_t point_t; + enum { no_children = SPTreeNode::no_children }; + +private: + 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); + + bool isCorrect() const; + void getAllIndices(unsigned int* indices) const; + unsigned int getDepth() const; + void computeNonEdgeForces(unsigned int point_index, double theta, double neg_f[], double* sum_Q) const; + void computeEdgeForces(unsigned int* row_P, unsigned int* col_P, double* val_P, int N, double* pos_f) const; + void print() const; + +private: + void fill(unsigned int N); + // Disallow copy + SPTree(const SPTree&) = delete; +}; + #endif From ef51dc0e6faffa60e3f0b3f44ff44d07dedc677e Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Fri, 1 Jun 2018 03:07:47 -0700 Subject: [PATCH 04/22] Namespace the internal implementation details. --- tsne/bh_sne_src/sptree.cpp | 227 ++++++++++++++++++------------------- tsne/bh_sne_src/sptree.h | 10 +- 2 files changed, 119 insertions(+), 118 deletions(-) diff --git a/tsne/bh_sne_src/sptree.cpp b/tsne/bh_sne_src/sptree.cpp index 6fe8dd6..58bff0f 100644 --- a/tsne/bh_sne_src/sptree.cpp +++ b/tsne/bh_sne_src/sptree.cpp @@ -44,6 +44,8 @@ using std::move; using std::unique_ptr; using std::vector; +namespace _sptree_internal { + template Cell::Cell(typename Cell::point_t&& p) : corner(move(p)) {} @@ -68,52 +70,10 @@ bool Cell::containsPoint(const double* point, const typename Cell: return true; } -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) { - // Compute mean, width, and height of current map (boundaries of SPTree) - int nD = 0; - point_t mean_Y, min_Y, max_Y; - mean_Y.fill(0.0); - min_Y.fill(DBL_MAX); - max_Y.fill(-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; - } - - double dbl_N = static_cast(N); - for(int d = 0; d < NDims; d++) mean_Y[d] /= dbl_N; - - // Construct SPTree - point_t width; - for(int d = 0; d < NDims; d++) { - width[d] = max(max_Y[d] - mean_Y[d], mean_Y[d] - min_Y[d]) + 1e-5; - } - widths.emplace_back(move(width)); - node.setCorner(move(mean_Y)); - fill(N); -} - // Constructor for SPTreeNode. template SPTreeNode::SPTreeNode(typename SPTreeNode::point_t&& inp_corner) : - boundary(move(inp_corner)) { - size = 0; - cum_size = 0; - - center_of_mass.fill(0.0); -} + boundary(move(inp_corner)) {} // Update the corner position. template @@ -126,12 +86,6 @@ bool SPTreeNode::is_leaf() const { return NDims == 0 || children[0].get() == nullptr; } -// Insert a point into the SPTree -template -bool SPTree::insert(unsigned int new_index) { - return node.insert(new_index, data, &widths, 1); -} - template bool SPTreeNode::insert(unsigned int new_index, const double* data, vector* widths, typename vector::size_type depth) { @@ -226,23 +180,6 @@ void SPTreeNode::subdivide(const double* data, vector* widths, t size = 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()); -} - - template bool SPTreeNode::isCorrect(const double* data, typename vector::const_iterator width) const { @@ -261,16 +198,6 @@ bool SPTreeNode::isCorrect(const double* data, return true; } - - -// Build a list of all indices in SPTree -template -void SPTree::getAllIndices(unsigned int* indices) const -{ - node.getAllIndices(indices, 0); -} - - // Build a list of all indices in SPTree template unsigned int SPTreeNode::getAllIndices(unsigned int* indices, unsigned int loc) const @@ -288,11 +215,6 @@ unsigned int SPTreeNode::getAllIndices(unsigned int* indices, unsigned in return loc; } -template -unsigned int SPTree::getDepth() const { - return node.getDepth(); -} - template unsigned int SPTreeNode::getDepth() const { if(is_leaf()) return 1; @@ -301,13 +223,6 @@ unsigned int SPTreeNode::getDepth() const { 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) const { - double max_width = maxWidth(); - node.computeNonEdgeForces(point_index, data + point_index * NDims, theta, neg_f, sum_Q, max_width*max_width); -} - // Compute non-edge forces using Barnes-Hut algorithm template void SPTreeNode::computeNonEdgeForces(unsigned int point_index, @@ -353,6 +268,115 @@ void SPTreeNode::computeNonEdgeForces(unsigned int point_index, } } +// 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 < 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 < 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(int i = 0; i < no_children; i++) children[i]->print(data); + } +} + +} // namespace _sptree_internal + +using namespace _sptree_internal; + +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) { + // Compute mean, width, and height of current map (boundaries of SPTree) + int nD = 0; + point_t mean_Y, min_Y, max_Y; + mean_Y.fill(0.0); + min_Y.fill(DBL_MAX); + max_Y.fill(-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; + } + + double dbl_N = static_cast(N); + for(int d = 0; d < NDims; d++) mean_Y[d] /= dbl_N; + + // Construct SPTree + point_t width; + for(int d = 0; d < NDims; d++) { + width[d] = max(max_Y[d] - mean_Y[d], mean_Y[d] - min_Y[d]) + 1e-5; + } + widths.emplace_back(move(width)); + node.setCorner(move(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, 1); +} + + +// 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(); + node.computeNonEdgeForces(point_index, data + point_index * NDims, theta, neg_f, sum_Q, max_width*max_width); +} // Computes edge forces template @@ -392,33 +416,6 @@ void SPTree::print() const { node.print(data); } -// 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 < 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 < 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(int i = 0; i < no_children; i++) children[i]->print(data); - } -} - // declare templates explicitly template class SPTree<2>; template class SPTree<3>; diff --git a/tsne/bh_sne_src/sptree.h b/tsne/bh_sne_src/sptree.h index 3e7a459..387af3a 100644 --- a/tsne/bh_sne_src/sptree.h +++ b/tsne/bh_sne_src/sptree.h @@ -38,6 +38,8 @@ #include #include +namespace _sptree_internal { + template class Cell final { public: @@ -116,15 +118,17 @@ struct SPTreeNode<0> enum { no_children = 1 }; }; +} // namespace _sptree_internal + template class SPTree { public: - typedef typename SPTreeNode::point_t point_t; - enum { no_children = SPTreeNode::no_children }; + typedef typename _sptree_internal::SPTreeNode::point_t point_t; + enum { no_children = _sptree_internal::SPTreeNode::no_children }; private: - SPTreeNode node; + _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. From 88c96ad4e55f032a402b70b8d0f358d270bc4959 Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Fri, 1 Jun 2018 04:07:47 -0700 Subject: [PATCH 05/22] Propagate constness in some arguments. Reduce manual memory management. Allow the compiler to take advantage of compile-time knowledge of the dimension when using the templated form of SPTree. Especially since, for us, D=2 or 3, the benefits of loop unrolling can be significant. Remove blas from the python build. It isn't actually being used, so there's no reason to introduce a dependency on it. Adjust the build settings for the python setup build to be the same as the Makefile, at least on linux. Specifically, turn on LTO and gc-sections. These make a large difference in binary size, especially when statically linking libc++. --- README.md | 1 - setup.py | 16 ++++++------- tsne/bh_sne_src/tsne.cpp | 51 ++++++++++++++++++---------------------- tsne/bh_sne_src/tsne.h | 12 +++++----- 4 files changed, 37 insertions(+), 43 deletions(-) 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..4352e4d 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: @@ -45,17 +45,17 @@ 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'], + extra_compile_args=['-msse2', '-O3', '-fPIC', '-w', '-ffast-math', '-std=c++11', + '-ffunction-sections', '-flto'], + extra_link_args=['-Wl,--gc-sections', '-flto'], 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'], + extra_compile_args=['-msse2', '-O3', '-fPIC', '-w', '-ffast-math', '-DTSNE3D', + '-std=c++11', '-ffunction-sections', '-flto'], + extra_link_args=['-Wl,--gc-sections', '-flto'], language='c++')] ext_modules = cythonize(ext_modules) diff --git a/tsne/bh_sne_src/tsne.cpp b/tsne/bh_sne_src/tsne.cpp index 5fac463..406b585 100644 --- a/tsne/bh_sne_src/tsne.cpp +++ b/tsne/bh_sne_src/tsne.cpp @@ -217,31 +217,29 @@ void TSNE::run(double* X, int N, int D, double* Y, int no_dims, double perplexit // 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) +void TSNE::computeGradient(double* P, unsigned int* inp_row_P, unsigned int* inp_col_P, double* inp_val_P, const double* Y, int N, int D, double* dC, double theta) { + assert(D == NDIMS); // 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 * NDIMS; + 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 * NDIMS, &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) { +void TSNE::computeExactGradient(double* P, const 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; @@ -291,7 +289,7 @@ void TSNE::computeExactGradient(double* P, double* Y, int N, int D, double* dC) // Evaluate t-SNE cost function (exactly) -double TSNE::evaluateError(double* P, double* Y, int N, int D) { +double TSNE::evaluateError(double* P, const double* Y, int N, int D) { // Compute the squared Euclidean distance matrix double* DD = (double*) malloc(N * N * sizeof(double)); @@ -327,40 +325,37 @@ double TSNE::evaluateError(double* P, double* Y, int N, int D) { } // 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) +double TSNE::evaluateError(unsigned int* row_P, unsigned int* col_P, double* val_P, const double* Y, int N, int D, double theta) { + assert(D == NDIMS); // Get estimate of normalization term - SPTree* tree = new SPTree(Y, N); - double* buff = (double*) calloc(D, sizeof(double)); + SPTree tree(Y, N); + double buff[NDIMS]; 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++) { - ind1 = n * D; + ind1 = n * NDIMS; 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]; + ind2 = col_P[i] * NDIMS; + for(int d = 0; d < NDIMS; d++) buff[d] = Y[ind1 + d]; + for(int d = 0; d < NDIMS; d++) buff[d] -= Y[ind2 + d]; + for(int d = 0; d < NDIMS; 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) { +void TSNE::computeGaussianPerplexity(const double* X, int N, int D, double* P, double perplexity) { // Compute the squared Euclidean distance matrix double* DD = (double*) malloc(N * N * sizeof(double)); @@ -618,7 +613,7 @@ void TSNE::symmetrizeMatrix(unsigned int** _row_P, unsigned int** _col_P, double } // Compute squared Euclidean distance matrix -void TSNE::computeSquaredEuclideanDistance(double* X, int N, int D, double* DD) { +void TSNE::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; diff --git a/tsne/bh_sne_src/tsne.h b/tsne/bh_sne_src/tsne.h index 93e203e..b65ff8a 100644 --- a/tsne/bh_sne_src/tsne.h +++ b/tsne/bh_sne_src/tsne.h @@ -50,14 +50,14 @@ class TSNE 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 computeGradient(double* P, unsigned int* inp_row_P, unsigned int* inp_col_P, double* inp_val_P, const double* Y, int N, int D, double* dC, double theta); + void computeExactGradient(double* P, const double* Y, int N, int D, double* dC); + double evaluateError(double* P, const double* Y, int N, int D); + double evaluateError(unsigned int* row_P, unsigned int* col_P, double* val_P, const 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(const 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); + void computeSquaredEuclideanDistance(const double* X, int N, int D, double* DD); double randn(); }; From 47eaa06ae008de53eb53cd10578ff8f455174ee6 Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Fri, 1 Jun 2018 04:30:18 -0700 Subject: [PATCH 06/22] More sorting and correction of includes. Also get rid of 'using namespace std' in tsne.cpp. Instead, 'using std::vector' gets us the same benefit with less risk of being broken by future versions of c++ introducing symbols with names we already use. Use nullptr instead of NULL. --- tsne/bh_sne_src/tsne.cpp | 24 ++++++++++++------------ tsne/bh_sne_src/vptree.h | 28 +++++++++++++--------------- 2 files changed, 25 insertions(+), 27 deletions(-) diff --git a/tsne/bh_sne_src/tsne.cpp b/tsne/bh_sne_src/tsne.cpp index 406b585..7a33f32 100644 --- a/tsne/bh_sne_src/tsne.cpp +++ b/tsne/bh_sne_src/tsne.cpp @@ -30,21 +30,21 @@ * */ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include "tsne.h" + #include "vptree.h" #include "sptree.h" -#include "tsne.h" -using namespace std; +#include +#include +#include +#include +#include +#include +#include +#include + +using std::vector; #ifdef TSNE3D #define NDIMS 3 diff --git a/tsne/bh_sne_src/vptree.h b/tsne/bh_sne_src/vptree.h index 9e301d6..ed64b2e 100644 --- a/tsne/bh_sne_src/vptree.h +++ b/tsne/bh_sne_src/vptree.h @@ -33,18 +33,16 @@ /* 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 +#include // for reverse, nth_element +#include // for DBL_MAX +#include // for sqrt +#include // for malloc, free, rand, RAND_MAX +#include // for priority_queue +#include // for vector + class DataPoint { int _ind; @@ -55,7 +53,7 @@ class DataPoint DataPoint() { _D = 1; _ind = -1; - _x = NULL; + _x = nullptr; } DataPoint(int D, int ind, double* x) { _D = D; @@ -71,10 +69,10 @@ class DataPoint for(int d = 0; d < _D; d++) _x[d] = other.x(d); } } - ~DataPoint() { if(_x != NULL) free(_x); } + ~DataPoint() { if(_x != nullptr) free(_x); } DataPoint& operator= (const DataPoint& other) { // asignment should free old object if(this != &other) { - if(_x != NULL) free(_x); + if(_x != nullptr) free(_x); _D = other.dimensionality(); _ind = other.index(); _x = (double*) malloc(_D * sizeof(double)); @@ -193,7 +191,7 @@ class VpTree Node* buildFromPoints( int lower, int upper ) { if (upper == lower) { // indicates that we're done here! - return NULL; + return nullptr; } // Lower index is center of current node @@ -229,7 +227,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,7 +240,7 @@ class VpTree } // Return if we arrived at a leaf - if(node->left == NULL && node->right == NULL) { + if(node->left == nullptr && node->right == nullptr) { return; } From 82917dd982272d218c58bd6f57c78042d727647c Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Mon, 4 Jun 2018 20:48:52 -0700 Subject: [PATCH 07/22] Fix some array init stuff. Add some assertions. --- tsne/bh_sne.pyx | 9 +++-- tsne/bh_sne_3d.pyx | 9 +++-- tsne/bh_sne_src/sptree.cpp | 80 ++++++++++++++++++++++---------------- tsne/bh_sne_src/sptree.h | 11 +++--- tsne/bh_sne_src/tsne.cpp | 20 +++------- tsne/bh_sne_src/tsne.h | 2 +- tsne/bh_sne_src/vptree.h | 44 ++++----------------- 7 files changed, 80 insertions(+), 95 deletions(-) diff --git a/tsne/bh_sne.pyx b/tsne/bh_sne.pyx index e45dd13..9a8b94b 100644 --- a/tsne/bh_sne.pyx +++ b/tsne/bh_sne.pyx @@ -21,8 +21,11 @@ cdef class BH_SNE: @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) + 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(d, 2) + assert(d, X.shape[0]) + assert(N, X.shape[1]) 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_3d.pyx b/tsne/bh_sne_3d.pyx index 315e107..1d19926 100644 --- a/tsne/bh_sne_3d.pyx +++ b/tsne/bh_sne_3d.pyx @@ -21,8 +21,11 @@ cdef class BH_SNE_3D: @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) + 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(d, 3) + assert(d, X.shape[0]) + assert(N, X.shape[1]) 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/sptree.cpp b/tsne/bh_sne_src/sptree.cpp index 58bff0f..780abef 100644 --- a/tsne/bh_sne_src/sptree.cpp +++ b/tsne/bh_sne_src/sptree.cpp @@ -36,18 +36,16 @@ #include // for DBL_MAX #include // for fprintf, stderr, size_t #include // for unique_ptr -#include // for move using std::max; using std::max_element; -using std::move; using std::unique_ptr; using std::vector; namespace _sptree_internal { template -Cell::Cell(typename Cell::point_t&& p) : corner(move(p)) {} +Cell::Cell(const typename Cell::point_t& p) : corner(p) {} template double Cell::getCorner(unsigned int d) const { @@ -55,15 +53,20 @@ double Cell::getCorner(unsigned int d) const { } template -void Cell::setCorner(typename Cell::point_t&& val) { +void Cell::setCorner(const typename Cell::point_t& val) { corner = val; } +template +void Cell::setCorner(unsigned int d, double val) { + corner[d] = val; +} + // Checks whether a point lies in a cell template bool Cell::containsPoint(const double* point, const typename Cell::point_t& width) const { - for(int d = 0; d < NDims; d++) { + 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; } @@ -72,13 +75,21 @@ bool Cell::containsPoint(const double* point, const typename Cell: // Constructor for SPTreeNode. template -SPTreeNode::SPTreeNode(typename SPTreeNode::point_t&& inp_corner) : - boundary(move(inp_corner)) {} +SPTreeNode::SPTreeNode(const typename SPTreeNode::point_t& inp_corner) : + boundary(inp_corner), size(0), cum_size(0) { + center_of_mass.fill(0.0); +} +// Constructor for SPTreeNode. +template +SPTreeNode::SPTreeNode() : + size(0), cum_size(0) { + center_of_mass.fill(0.0); +} // Update the corner position. template -void SPTreeNode::setCorner(typename SPTreeNode::point_t&& inp_corner) { - boundary.setCorner(move(inp_corner)); +void SPTreeNode::setCorner(const typename SPTreeNode::point_t& inp_corner) { + boundary.setCorner(inp_corner); } template @@ -91,8 +102,7 @@ bool SPTreeNode::insert(unsigned int new_index, const double* data, vecto { // Ignore objects which do not belong in this quad tree const double* point = data + new_index * NDims; - const auto& width = (*widths)[depth]; - if(!boundary.containsPoint(point, width)) + if(!boundary.containsPoint(point, (*widths)[depth])) return false; // Online update of cumulative size and center-of-mass @@ -144,27 +154,27 @@ void SPTreeNode::subdivide(const double* data, vector* widths, t // If nessessary, add to the width. if (depth+1 == widths->size()) { // extend the list. - point_t child_widths; + 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]; } - widths->emplace_back(move(child_widths)); } const point_t& new_width = (*widths)[depth+1]; // Create new children for(unsigned int i = 0; i < no_children; i++) { unsigned int div = 1; - point_t new_corner; + children[i].reset(new SPTreeNode()); + Cell& new_corner = children[i]->boundary; for(unsigned int d = 0; d < NDims; d++) { if((i / div) % 2 == 1) - new_corner[d] = boundary.getCorner(d) - new_width[d]; + new_corner.setCorner(d, boundary.getCorner(d) - new_width[d]); else - new_corner[d] = boundary.getCorner(d) + new_width[d]; + new_corner.setCorner(d, boundary.getCorner(d) + new_width[d]); div *= 2; } - children[i].reset(new SPTreeNode(move(new_corner))); } // Move existing points to correct children @@ -231,7 +241,10 @@ void SPTreeNode::computeNonEdgeForces(unsigned int point_index, { // 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 || + (size == 1 && + __builtin_expect(index[0] == point_index, 0) && + is_leaf())) return; // Compute distance between point and center-of-mass double sqdist = .0; @@ -306,7 +319,7 @@ double SPTree::maxWidth() const { // Top-node constructor for SPTree -- build tree, too! template -SPTree::SPTree(const double* inp_data, unsigned int N) { +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; point_t mean_Y, min_Y, max_Y; @@ -314,32 +327,33 @@ SPTree::SPTree(const double* inp_data, unsigned int N) { min_Y.fill(DBL_MAX); max_Y.fill(-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]; + const double *elem = inp_data; + for(unsigned int n = 0; n < N; ++n) { + for(unsigned 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]; } - nD += NDims; + elem += NDims; } double dbl_N = static_cast(N); for(int d = 0; d < NDims; d++) mean_Y[d] /= dbl_N; // Construct SPTree - point_t width; + 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; } - widths.emplace_back(move(width)); - node.setCorner(move(mean_Y)); + 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, 1); + return node.insert(new_index, data, &widths, 0); } @@ -393,17 +407,17 @@ void SPTree::computeEdgeForces(unsigned int* row_P, unsigned int* col_P, sqdist = 1.0; unsigned int ind2 = col_P[i] * NDims; + std::array diffs; for(unsigned int d = 0; d < NDims; d++) { - double diff = data[ind1 + d] - data[ind2 + d]; - sqdist += diff * diff; + diffs[d] = data[ind1 + d] - data[ind2 + d]; + sqdist += diffs[d] * diffs[d]; } sqdist = val_P[i] / sqdist; // Sum positive force for(unsigned int d = 0; d < NDims; d++) { - double diff = data[ind1 + d] - data[ind2 + d]; - pos_f[ind1 + d] += sqdist * diff; + pos_f[ind1 + d] += sqdist * diffs[d]; } } ind1 += NDims; diff --git a/tsne/bh_sne_src/sptree.h b/tsne/bh_sne_src/sptree.h index 387af3a..94000ca 100644 --- a/tsne/bh_sne_src/sptree.h +++ b/tsne/bh_sne_src/sptree.h @@ -46,10 +46,11 @@ class Cell final { typedef std::array point_t; Cell() = default; - Cell(point_t&&); + explicit Cell(const point_t&); double getCorner(unsigned int d) const; - void setCorner(point_t&& inp_corner); + void setCorner(const point_t& inp_corner); + void setCorner(unsigned int d, double val); bool containsPoint(const double* point, const point_t& width) const; private: @@ -93,10 +94,10 @@ class SPTreeNode final { bool is_leaf() const; public: - SPTreeNode() = default; - SPTreeNode(point_t&& corner); + SPTreeNode(); + explicit SPTreeNode(const point_t& corner); - void setCorner(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); bool isCorrect(const double* data, typename std::vector::const_iterator width) const; diff --git a/tsne/bh_sne_src/tsne.cpp b/tsne/bh_sne_src/tsne.cpp index 7a33f32..51ed577 100644 --- a/tsne/bh_sne_src/tsne.cpp +++ b/tsne/bh_sne_src/tsne.cpp @@ -83,12 +83,9 @@ 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"); @@ -164,8 +161,8 @@ void TSNE::run(double* X, int N, int D, double* Y, int no_dims, double perplexit 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, Y, N, no_dims, dY.data()); + else computeGradient(P, row_P, col_P, val_P, Y, N, no_dims, 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); @@ -203,9 +200,6 @@ void TSNE::run(double* X, int N, int D, double* Y, int no_dims, double perplexit 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; @@ -219,8 +213,6 @@ void TSNE::run(double* X, int N, int D, double* Y, int no_dims, double perplexit // 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, const double* Y, int N, int D, double* dC, double theta) { - assert(D == NDIMS); - // Construct space-partitioning tree on current map SPTree tree(Y, N); @@ -426,7 +418,7 @@ void TSNE::computeGaussianPerplexity(const double* X, int N, int D, double* P, d // 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) { +void TSNE::computeGaussianPerplexity(const double* X, int N, int D, unsigned int** _row_P, unsigned int** _col_P, double** _val_P, double perplexity, int K) { if(perplexity > K) fprintf(stderr,"Perplexity should be lower than K!\n"); diff --git a/tsne/bh_sne_src/tsne.h b/tsne/bh_sne_src/tsne.h index b65ff8a..200402b 100644 --- a/tsne/bh_sne_src/tsne.h +++ b/tsne/bh_sne_src/tsne.h @@ -56,7 +56,7 @@ class TSNE double evaluateError(unsigned int* row_P, unsigned int* col_P, double* val_P, const double* Y, int N, int D, double theta); void zeroMean(double* X, int N, int D); void computeGaussianPerplexity(const 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 computeGaussianPerplexity(const double* X, int N, int D, unsigned int** _row_P, unsigned int** _col_P, double** _val_P, double perplexity, int K); void computeSquaredEuclideanDistance(const double* X, int N, int D, double* DD); double randn(); }; diff --git a/tsne/bh_sne_src/vptree.h b/tsne/bh_sne_src/vptree.h index ed64b2e..a8bd7d7 100644 --- a/tsne/bh_sne_src/vptree.h +++ b/tsne/bh_sne_src/vptree.h @@ -48,49 +48,21 @@ class DataPoint int _ind; public: - double* _x; - int _D; - DataPoint() { - _D = 1; - _ind = -1; - _x = nullptr; - } - 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 != nullptr) free(_x); } - DataPoint& operator= (const DataPoint& other) { // asignment should free old object - if(this != &other) { - if(_x != nullptr) 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; - } + std::vector _x; + DataPoint() : _ind(-1) {} + DataPoint(int D, int ind, const double* x) : _ind(ind), _x(x, x+D) {} + DataPoint(const DataPoint& other) = default; int index() const { return _ind; } - int dimensionality() const { return _D; } + int dimensionality() const { return (int)_x.size(); } 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; + const std::vector& x1 = t1._x; + const std::vector& x2 = t2._x; double diff; - for(int d = 0; d < t1._D; d++) { + for(int d = 0; d < x1.size(); d++) { diff = (x1[d] - x2[d]); dd += diff * diff; } From dc8199c65840248be1944846a6cf7fe8c35fe598 Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Wed, 13 Jun 2018 17:09:35 -0700 Subject: [PATCH 08/22] Factor out main. Because we use tsne through the cython wrapper, we don't actually use the c++ main method. --- tsne/bh_sne_src/Makefile | 7 +++- tsne/bh_sne_src/main.cpp | 81 ++++++++++++++++++++++++++++++++++++++++ tsne/bh_sne_src/tsne.cpp | 48 ------------------------ 3 files changed, 86 insertions(+), 50 deletions(-) create mode 100644 tsne/bh_sne_src/main.cpp diff --git a/tsne/bh_sne_src/Makefile b/tsne/bh_sne_src/Makefile index 20b37f2..d6bebef 100644 --- a/tsne/bh_sne_src/Makefile +++ b/tsne/bh_sne_src/Makefile @@ -10,10 +10,10 @@ LDFLAGS += -Wl,--gc-sections all: bh_tsne bh_tsne_3d -bh_tsne: tsne.o sptree.o +bh_tsne: main.o tsne.o sptree.o $(CXX) $(CFLAGS) $(LDFLAGS) $^ -o $@ -bh_tsne_3d: tsne_3d.o sptree.o +bh_tsne_3d: tsne_3d.o sptree.o main.o $(CXX) $(CFLAGS) $(LDFLAGS) $^ -o $@ sptree.o: sptree.cpp sptree.h @@ -25,5 +25,8 @@ tsne.o: tsne.cpp tsne.h sptree.h vptree.h tsne_3d.o: tsne.cpp tsne.h sptree.h vptree.h $(CXX) $(CFLAGS) -DTSNE3D -c $< -o $@ +main.o: main.cpp tsne.h sptree.h vptree.h + $(CXX) $(CFLAGS) -c $< -o $@ + clean: rm -Rf *.o bh_tsne bh_tsne_3d diff --git a/tsne/bh_sne_src/main.cpp b/tsne/bh_sne_src/main.cpp new file mode 100644 index 0000000..e7e159b --- /dev/null +++ b/tsne/bh_sne_src/main.cpp @@ -0,0 +1,81 @@ +/* + * + * 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; + TSNE tsne; + + // 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/tsne.cpp b/tsne/bh_sne_src/tsne.cpp index 51ed577..2cb7bea 100644 --- a/tsne/bh_sne_src/tsne.cpp +++ b/tsne/bh_sne_src/tsne.cpp @@ -727,51 +727,3 @@ void TSNE::save_csv(const char* csv_file, double* Y, int N, int D) { csv.close(); } - -// Function that runs the Barnes-Hut implementation of t-SNE -int main(int argc, char *argv[]) { - - // 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]; - } - - 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; - } - delete(tsne); -} From 0019c9baea11ac9bfec47f5e2853ee23a8bd68e8 Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Wed, 13 Jun 2018 18:08:04 -0700 Subject: [PATCH 09/22] Use dynamic dispatch for choosing no_dims. As opposed to tricks with linking and #defines which cython doesn't do a great job obeying. --- tsne/bh_sne_src/Makefile | 12 +-- tsne/bh_sne_src/sptree.cpp | 12 +-- tsne/bh_sne_src/tsne.cpp | 202 ++++++++++++++++++++++++------------- tsne/bh_sne_src/tsne.h | 12 --- 4 files changed, 142 insertions(+), 96 deletions(-) diff --git a/tsne/bh_sne_src/Makefile b/tsne/bh_sne_src/Makefile index d6bebef..85c8902 100644 --- a/tsne/bh_sne_src/Makefile +++ b/tsne/bh_sne_src/Makefile @@ -8,25 +8,19 @@ CFLAGS += -std=c++11 -ffast-math -O3 -flto CFLAGS += -ffunction-sections LDFLAGS += -Wl,--gc-sections -all: bh_tsne bh_tsne_3d +all: bh_tsne bh_tsne: main.o tsne.o sptree.o $(CXX) $(CFLAGS) $(LDFLAGS) $^ -o $@ -bh_tsne_3d: tsne_3d.o sptree.o main.o - $(CXX) $(CFLAGS) $(LDFLAGS) $^ -o $@ - sptree.o: sptree.cpp sptree.h $(CXX) $(CFLAGS) -c $< -o $@ -tsne.o: tsne.cpp tsne.h sptree.h vptree.h +main.o tsne.o: tsne.cpp tsne.h sptree.h vptree.h $(CXX) $(CFLAGS) -c $< -o $@ -tsne_3d.o: tsne.cpp tsne.h sptree.h vptree.h - $(CXX) $(CFLAGS) -DTSNE3D -c $< -o $@ - main.o: main.cpp tsne.h sptree.h vptree.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/sptree.cpp b/tsne/bh_sne_src/sptree.cpp index 780abef..95f8811 100644 --- a/tsne/bh_sne_src/sptree.cpp +++ b/tsne/bh_sne_src/sptree.cpp @@ -37,6 +37,7 @@ #include // for fprintf, stderr, size_t #include // for unique_ptr +using std::array; using std::max; using std::max_element; using std::unique_ptr; @@ -249,9 +250,10 @@ void SPTreeNode::computeNonEdgeForces(unsigned int point_index, // Compute distance between point and center-of-mass double sqdist = .0; + array diff; for(int i = 0; i < NDims; ++i) { - double diff = data_point[i] - center_of_mass[i]; - sqdist += diff * diff; + diff[i] = data_point[i] - center_of_mass[i]; + sqdist += diff[i] * diff[i]; } // Check whether we can use this node as a "summary" @@ -264,11 +266,7 @@ void SPTreeNode::computeNonEdgeForces(unsigned int point_index, *sum_Q += mult; mult *= sqdist; for(size_t d = 0; d < NDims; ++d) { - // recompute here rather than storing from before because memory - // locality matters more than an extra couple of additions and a - // subtraction. - double diff = data_point[d] - center_of_mass[d]; - neg_f[d] += mult * diff; + neg_f[d] += mult * diff[d]; } } else { diff --git a/tsne/bh_sne_src/tsne.cpp b/tsne/bh_sne_src/tsne.cpp index 2cb7bea..0110518 100644 --- a/tsne/bh_sne_src/tsne.cpp +++ b/tsne/bh_sne_src/tsne.cpp @@ -44,16 +44,33 @@ #include #include +using std::array; using std::vector; -#ifdef TSNE3D -#define NDIMS 3 -#else -#define NDIMS 2 -#endif +namespace { + +void symmetrizeMatrix(unsigned int** row_P, unsigned int** col_P, double** 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, unsigned int** _row_P, unsigned int** _col_P, double** _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 -void TSNE::run(double* X, int N, int D, double* Y, int no_dims, double perplexity, double theta, int rand_seed, +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 ) { @@ -161,8 +178,8 @@ void TSNE::run(double* X, int N, int D, double* Y, int no_dims, double perplexit for(int iter = 0; iter < max_iter; iter++) { // Compute (approximate) gradient - if(exact) computeExactGradient(P, Y, N, no_dims, dY.data()); - else computeGradient(P, row_P, col_P, val_P, Y, N, no_dims, dY.data(), theta); + if(exact) computeExactGradient(P, Y, N, dY.data()); + else computeGradient(P, row_P, col_P, val_P, 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); @@ -173,7 +190,7 @@ void TSNE::run(double* X, int N, int D, double* Y, int no_dims, double perplexit 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) { @@ -186,8 +203,8 @@ void TSNE::run(double* X, int N, int D, double* Y, int no_dims, double perplexit 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(exact) C = evaluateError(P, Y, N); + else C = evaluateError(row_P, col_P, val_P, Y, N, theta); // doing approximate computation here! if(iter == 0) fprintf(stderr,"Iteration %d: error is %f\n", iter + 1, C); else { @@ -211,18 +228,19 @@ void TSNE::run(double* X, int N, int D, double* Y, int no_dims, double perplexit // 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, const 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(Y, N); + SPTree tree(Y, N); // Compute all terms required for t-SNE gradient double sum_Q = .0; - auto len = N * NDIMS; + 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 * NDIMS, &sum_Q); + 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 < len; i++) { @@ -231,19 +249,18 @@ void TSNE::computeGradient(double* P, unsigned int* inp_row_P, unsigned int* inp } // Compute gradient of the t-SNE cost function (exact) -void TSNE::computeExactGradient(double* P, const double* Y, int N, int D, double* dC) { +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++) { @@ -273,25 +290,21 @@ void TSNE::computeExactGradient(double* P, const double* Y, int N, int D, double 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, const 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; + vector Q(N * N); for(int n = 0; n < N; n++) { for(int m = 0; m < N; m++) { if(n != m) { @@ -310,20 +323,16 @@ double TSNE::evaluateError(double* P, const double* Y, int N, int D) { C += P[n] * log((P[n] + FLT_MIN) / (Q[n] + FLT_MIN)); } - // Clean up memory - free(DD); - free(Q); return C; } // Evaluate t-SNE cost function (approximately) -double TSNE::evaluateError(unsigned int* row_P, unsigned int* col_P, double* val_P, const 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) { - assert(D == NDIMS); - // Get estimate of normalization term - SPTree tree(Y, N); - double buff[NDIMS]; + 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); @@ -331,13 +340,13 @@ double TSNE::evaluateError(unsigned int* row_P, unsigned int* col_P, double* val int ind1, ind2; double C = .0, Q; for(int n = 0; n < N; n++) { - ind1 = n * NDIMS; + ind1 = n * D; for(int i = row_P[n]; i < row_P[n + 1]; i++) { Q = .0; - ind2 = col_P[i] * NDIMS; - for(int d = 0; d < NDIMS; d++) buff[d] = Y[ind1 + d]; - for(int d = 0; d < NDIMS; d++) buff[d] -= Y[ind2 + d]; - for(int d = 0; d < NDIMS; d++) Q += buff[d] * buff[d]; + 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]; Q = (1.0 / (1.0 + Q)) / sum_Q; C += val_P[i] * log((val_P[i] + FLT_MIN) / (Q + FLT_MIN)); } @@ -347,12 +356,11 @@ double TSNE::evaluateError(unsigned int* row_P, unsigned int* col_P, double* val // Compute input similarities with a fixed perplexity -void TSNE::computeGaussianPerplexity(const double* X, int N, int D, double* P, double perplexity) { +void computeGaussianPerplexity(const 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); + vector DD(N * N); + computeSquaredEuclideanDistance(X, N, D, DD.data()); // Compute the Gaussian kernel row by row int nN = 0; @@ -411,14 +419,11 @@ void TSNE::computeGaussianPerplexity(const double* X, int N, int D, double* P, d for(int m = 0; m < N; m++) P[nN + m] /= sum_P; nN += N; } - - // Clean up memory - free(DD); DD = NULL; } // Compute input similarities with a fixed perplexity using ball trees (this function allocates memory another function should free) -void TSNE::computeGaussianPerplexity(const double* X, int N, int D, unsigned int** _row_P, unsigned int** _col_P, double** _val_P, double perplexity, int K) { +void computeGaussianPerplexity(const double* X, int N, int D, unsigned int** _row_P, unsigned int** _col_P, double** _val_P, double perplexity, int K) { if(perplexity > K) fprintf(stderr,"Perplexity should be lower than K!\n"); @@ -430,16 +435,15 @@ void TSNE::computeGaussianPerplexity(const double* X, int N, int D, unsigned int 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); } + 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; // Build ball tree on data set - VpTree* tree = new VpTree(); + VpTree tree; 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); + tree.create(obj_X); // Loop over all points to find nearest neighbors fprintf(stderr,"Building tree...\n"); @@ -452,7 +456,7 @@ void TSNE::computeGaussianPerplexity(const double* X, int N, int D, unsigned int // 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; @@ -508,16 +512,11 @@ void TSNE::computeGaussianPerplexity(const double* X, int N, int D, unsigned int 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(unsigned int** _row_P, unsigned int** _col_P, double** _val_P, int N) { // Get sparse matrix unsigned int* row_P = *_row_P; @@ -605,7 +604,25 @@ void TSNE::symmetrizeMatrix(unsigned int** _row_P, unsigned int** _col_P, double } // Compute squared Euclidean distance matrix -void TSNE::computeSquaredEuclideanDistance(const 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) { + 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; + } + } +} + +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; @@ -624,11 +641,38 @@ void TSNE::computeSquaredEuclideanDistance(const double* X, int N, int D, double // Makes data zero-mean -void TSNE::zeroMean(double* X, int N, int D) { +void zeroMean(double* X, int N, int D) { + + // Compute data mean + vector mean(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++) { + mean[d] /= (double) N; + } + + // Subtract data mean + nD = 0; + for(int n = 0; n < N; n++) { + for(int d = 0; d < D; d++) { + X[nD + d] -= mean[d]; + } + nD += D; + } +} + +// Makes data zero-mean +template +void zeroMean(double* X, int N) { // Compute data mean - double* mean = (double*) calloc(D, sizeof(double)); - if(mean == NULL) { fprintf(stderr,"Memory allocation failed!\n"); exit(1); } + array mean; + mean.fill(0); int nD = 0; for(int n = 0; n < N; n++) { for(int d = 0; d < D; d++) { @@ -648,12 +692,11 @@ void TSNE::zeroMean(double* X, int N, int D) { } nD += D; } - free(mean); mean = NULL; } // Generates a Gaussian random number -double TSNE::randn() { +double randn() { double x, y, radius; do { x = 2 * (rand() / ((double) RAND_MAX + 1)) - 1; @@ -666,6 +709,29 @@ double TSNE::randn() { return x; } +} // 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 + ) { + 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"); + } +} + // 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) { diff --git a/tsne/bh_sne_src/tsne.h b/tsne/bh_sne_src/tsne.h index 200402b..2b2597e 100644 --- a/tsne/bh_sne_src/tsne.h +++ b/tsne/bh_sne_src/tsne.h @@ -46,19 +46,7 @@ class TSNE ); 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); - -private: - void computeGradient(double* P, unsigned int* inp_row_P, unsigned int* inp_col_P, double* inp_val_P, const double* Y, int N, int D, double* dC, double theta); - void computeExactGradient(double* P, const double* Y, int N, int D, double* dC); - double evaluateError(double* P, const double* Y, int N, int D); - double evaluateError(unsigned int* row_P, unsigned int* col_P, double* val_P, const double* Y, int N, int D, double theta); - void zeroMean(double* X, int N, int D); - void computeGaussianPerplexity(const double* X, int N, int D, double* P, double perplexity); - void computeGaussianPerplexity(const double* X, int N, int D, unsigned int** _row_P, unsigned int** _col_P, double** _val_P, double perplexity, int K); - void computeSquaredEuclideanDistance(const double* X, int N, int D, double* DD); - double randn(); }; #endif From 930525e2a6d09de387fcf176df0053a042887796 Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Wed, 13 Jun 2018 18:23:10 -0700 Subject: [PATCH 10/22] Use a namespace instead of class. The class had no members, which made it useless. --- tsne/bh_sne_src/main.cpp | 7 +++--- tsne/bh_sne_src/tsne.cpp | 52 +++++++++++++++++++--------------------- tsne/bh_sne_src/tsne.h | 25 +++++++++---------- 3 files changed, 41 insertions(+), 43 deletions(-) diff --git a/tsne/bh_sne_src/main.cpp b/tsne/bh_sne_src/main.cpp index e7e159b..86d8205 100644 --- a/tsne/bh_sne_src/main.cpp +++ b/tsne/bh_sne_src/main.cpp @@ -57,23 +57,22 @@ int main(int argc, char *argv[]) { double perc_landmarks; double perplexity, theta, *data; int rand_seed = -1; - TSNE tsne; // Read the parameters and the dataset - if(tsne.load_data(dat_file, &data, &origN, &D, &no_dims, &theta, &perplexity, &rand_seed, &max_iter)) { + 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); + 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); + 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/tsne.cpp b/tsne/bh_sne_src/tsne.cpp index 0110518..cdf1c95 100644 --- a/tsne/bh_sne_src/tsne.cpp +++ b/tsne/bh_sne_src/tsne.cpp @@ -47,8 +47,11 @@ using std::array; using std::vector; +namespace TSNE { namespace { +static inline double sign(double x) { return (x == .0 ? .0 : (x < .0 ? -1.0 : 1.0)); } + void symmetrizeMatrix(unsigned int** row_P, unsigned int** col_P, double** 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); @@ -115,14 +118,14 @@ void run(double* X, int N, int D, double* Y, double perplexity, double theta, in 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; + vector P; + unsigned int* row_P; unsigned int* col_P; double* 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); + computeGaussianPerplexity(X, N, D, P.data(), perplexity); // Symmetrize input similarities fprintf(stderr,"Symmetrizing...\n"); @@ -178,8 +181,8 @@ void run(double* X, int N, int D, double* Y, double perplexity, double theta, in for(int iter = 0; iter < max_iter; iter++) { // Compute (approximate) gradient - if(exact) computeExactGradient(P, Y, N, dY.data()); - else computeGradient(P, row_P, col_P, val_P, Y, N, dY.data(), theta); + if(exact) computeExactGradient(P.data(), Y, N, dY.data()); + else computeGradient(P.data(), row_P, col_P, val_P, 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); @@ -203,7 +206,7 @@ void run(double* X, int N, int D, double* Y, double perplexity, double theta, in if (iter > 0 && (iter % 50 == 0 || iter == max_iter - 1)) { end = clock(); double C = .0; - if(exact) C = evaluateError(P, Y, N); + if(exact) C = evaluateError(P.data(), Y, N); else C = evaluateError(row_P, col_P, val_P, Y, N, theta); // doing approximate computation here! if(iter == 0) fprintf(stderr,"Iteration %d: error is %f\n", iter + 1, C); @@ -217,8 +220,7 @@ void run(double* X, int N, int D, double* Y, double perplexity, double theta, in end = clock(); total_time += (float) (end - start) / CLOCKS_PER_SEC; // Clean up memory - if(exact) free(P); - else { + if(!exact) { free(row_P); row_P = NULL; free(col_P); col_P = NULL; free(val_P); val_P = NULL; @@ -524,8 +526,7 @@ void symmetrizeMatrix(unsigned int** _row_P, unsigned int** _col_P, double** _va double* 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); } + vector row_counts(N); for(int n = 0; n < N; n++) { for(int i = row_P[n]; i < row_P[n + 1]; i++) { @@ -555,8 +556,7 @@ void symmetrizeMatrix(unsigned int** _row_P, unsigned int** _col_P, double** _va 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); } + 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]) @@ -597,10 +597,6 @@ void symmetrizeMatrix(unsigned int** _row_P, unsigned int** _col_P, double** _va 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; } // Compute squared Euclidean distance matrix @@ -712,20 +708,20 @@ double randn() { } // 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, +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); + 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); + 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"); @@ -734,7 +730,7 @@ void TSNE::run(double* X, int N, int D, double* Y, int no_dims, double perplexit // 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) { +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; @@ -758,7 +754,7 @@ bool TSNE::load_data(const char* dat_file, double** data, int* n, int* d, int* n } // 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) { +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; @@ -775,7 +771,7 @@ void TSNE::save_data(const char* res_file, double* data, int* landmarks, double* 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) { +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++) { @@ -793,3 +789,5 @@ void TSNE::save_csv(const char* csv_file, double* Y, int N, int D) { csv.close(); } + +} // namespace TSNE diff --git a/tsne/bh_sne_src/tsne.h b/tsne/bh_sne_src/tsne.h index 2b2597e..caa4fbf 100644 --- a/tsne/bh_sne_src/tsne.h +++ b/tsne/bh_sne_src/tsne.h @@ -34,19 +34,20 @@ #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 +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); -}; +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); + +} // namespace TSNE #endif From 6885a1930abe2c96f628e80d0de751639122e070 Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Wed, 13 Jun 2018 18:37:46 -0700 Subject: [PATCH 11/22] Remove obsolete bh_sne_3d. --- setup.py | 8 -------- tsne/bh_sne.pyx | 18 +++--------------- tsne/bh_sne_3d.pyx | 31 ------------------------------- 3 files changed, 3 insertions(+), 54 deletions(-) delete mode 100644 tsne/bh_sne_3d.pyx diff --git a/setup.py b/setup.py index 4352e4d..ea3a463 100644 --- a/setup.py +++ b/setup.py @@ -48,14 +48,6 @@ extra_compile_args=['-msse2', '-O3', '-fPIC', '-w', '-ffast-math', '-std=c++11', '-ffunction-sections', '-flto'], extra_link_args=['-Wl,--gc-sections', '-flto'], - 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/'], - extra_compile_args=['-msse2', '-O3', '-fPIC', '-w', '-ffast-math', '-DTSNE3D', - '-std=c++11', '-ffunction-sections', '-flto'], - extra_link_args=['-Wl,--gc-sections', '-flto'], language='c++')] ext_modules = cythonize(ext_modules) diff --git a/tsne/bh_sne.pyx b/tsne/bh_sne.pyx index 9a8b94b..33cc7a5 100644 --- a/tsne/bh_sne.pyx +++ b/tsne/bh_sne.pyx @@ -4,28 +4,16 @@ 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 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: - 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, 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(d, 2) - assert(d, X.shape[0]) assert(N, X.shape[1]) - 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) + 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 1d19926..0000000 --- a/tsne/bh_sne_3d.pyx +++ /dev/null @@ -1,31 +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, 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(d, 3) - assert(d, X.shape[0]) - assert(N, X.shape[1]) - 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 From 6da8ab83c2e4532f80053c071aaa7fe240ee252a Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Wed, 20 Jun 2018 10:46:30 -0700 Subject: [PATCH 12/22] Break out io stuff into separate obj. It's not used by the python wrapper, so no need to compile it. --- tsne/bh_sne_src/Makefile | 7 ++- tsne/bh_sne_src/io.cpp | 103 +++++++++++++++++++++++++++++++++++++++ tsne/bh_sne_src/tsne.cpp | 70 ++------------------------ 3 files changed, 113 insertions(+), 67 deletions(-) create mode 100644 tsne/bh_sne_src/io.cpp diff --git a/tsne/bh_sne_src/Makefile b/tsne/bh_sne_src/Makefile index 85c8902..42b4bcf 100644 --- a/tsne/bh_sne_src/Makefile +++ b/tsne/bh_sne_src/Makefile @@ -10,17 +10,20 @@ LDFLAGS += -Wl,--gc-sections all: bh_tsne -bh_tsne: main.o tsne.o sptree.o +bh_tsne: main.o io.o tsne.o sptree.o $(CXX) $(CFLAGS) $(LDFLAGS) $^ -o $@ sptree.o: sptree.cpp sptree.h $(CXX) $(CFLAGS) -c $< -o $@ -main.o tsne.o: tsne.cpp tsne.h sptree.h vptree.h +tsne.o: tsne.cpp tsne.h sptree.h vptree.h $(CXX) $(CFLAGS) -c $< -o $@ 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 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/tsne.cpp b/tsne/bh_sne_src/tsne.cpp index cdf1c95..6a3b762 100644 --- a/tsne/bh_sne_src/tsne.cpp +++ b/tsne/bh_sne_src/tsne.cpp @@ -30,11 +30,10 @@ * */ -#include "tsne.h" -#include "vptree.h" -#include "sptree.h" +#include "tsne.h" +#include #include #include #include @@ -44,6 +43,9 @@ #include #include +#include "vptree.h" +#include "sptree.h" + using std::array; using std::vector; @@ -728,66 +730,4 @@ void run(double* X, int N, int D, double* Y, int no_dims, double perplexity, dou } } -// 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 From c2a2ef124452641cb359e1df28daba488a538855 Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Wed, 20 Jun 2018 18:51:54 -0700 Subject: [PATCH 13/22] Improve build/link. Properly take advantage of dynamic dispatch. --- Makefile | 8 +++++--- setup.py | 6 +++--- tsne/__init__.py | 10 +--------- tsne/bh_sne.pyx | 31 +++++++++++++++++++++++++------ 4 files changed, 34 insertions(+), 21 deletions(-) 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/setup.py b/setup.py index ea3a463..88b1447 100644 --- a/setup.py +++ b/setup.py @@ -44,10 +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/'], + include_dirs=[numpy.get_include(), 'tsne/bh_sne_src/'], extra_compile_args=['-msse2', '-O3', '-fPIC', '-w', '-ffast-math', '-std=c++11', - '-ffunction-sections', '-flto'], - extra_link_args=['-Wl,--gc-sections', '-flto'], + '-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 33cc7a5..77e4832 100644 --- a/tsne/bh_sne.pyx +++ b/tsne/bh_sne.pyx @@ -5,15 +5,34 @@ cimport cython from libcpp cimport bool cdef extern from "tsne.h" namespace "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) + 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: @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, 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') + @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]) - 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) + 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 From b266c94df60c955afbe83c3e90c1e0a21f72f597 Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Wed, 20 Jun 2018 18:52:48 -0700 Subject: [PATCH 14/22] Remove size from SPTreeNode. --- tsne/bh_sne_src/sptree.cpp | 165 ++++++++++++++++++++++--------------- tsne/bh_sne_src/sptree.h | 15 ++-- tsne/bh_sne_src/tsne.cpp | 17 ++-- 3 files changed, 120 insertions(+), 77 deletions(-) diff --git a/tsne/bh_sne_src/sptree.cpp b/tsne/bh_sne_src/sptree.cpp index 95f8811..10547de 100644 --- a/tsne/bh_sne_src/sptree.cpp +++ b/tsne/bh_sne_src/sptree.cpp @@ -43,6 +43,7 @@ using std::max_element; using std::unique_ptr; using std::vector; +namespace TSNE { namespace _sptree_internal { template @@ -77,13 +78,13 @@ bool Cell::containsPoint(const double* point, const typename Cell: // Constructor for SPTreeNode. template SPTreeNode::SPTreeNode(const typename SPTreeNode::point_t& inp_corner) : - boundary(inp_corner), size(0), cum_size(0) { + boundary(inp_corner), cum_size(0) { center_of_mass.fill(0.0); } // Constructor for SPTreeNode. template SPTreeNode::SPTreeNode() : - size(0), cum_size(0) { + cum_size(0) { center_of_mass.fill(0.0); } @@ -95,7 +96,20 @@ void SPTreeNode::setCorner(const typename SPTreeNode::point_t& inp template bool SPTreeNode::is_leaf() const { - return NDims == 0 || children[0].get() == nullptr; + return NDims == 0 || cum_size <= QT_NODE_CAPACITY; +} + +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; } template @@ -107,47 +121,48 @@ bool SPTreeNode::insert(unsigned int new_index, const double* data, vecto 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(size < QT_NODE_CAPACITY && is_leaf()) { + if(size < QT_NODE_CAPACITY) { index[size] = new_index; - size++; return true; + } 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; + for(unsigned int d = 0; d < NDims; d++) { + if(__builtin_expect(point[d] != data[index[n] * NDims + d], 1)) { duplicate = false; break; } + } + if (__builtin_expect(duplicate, 0)) { + cum_size--; + return true; + } + } // Otherwise, we need to subdivide the current cell + subdivide(data, widths, depth); } - // Don't add duplicates for now (this is not very nice) - 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; } - } - if (duplicate) { - return true; - } - } - - // Otherwise, we need to subdivide the current cell - if(is_leaf()) subdivide(data, widths, depth); - // Find out where the point can be inserted - for(auto& child : children) { - if(child->insert(new_index, data, widths, depth+1)) { - 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 SPTreeNode::subdivide(const double* data, vector* widths, typename vector::size_type depth) { @@ -166,39 +181,32 @@ void SPTreeNode::subdivide(const double* data, vector* widths, t // Create new children for(unsigned int i = 0; i < no_children; i++) { - unsigned int div = 1; children[i].reset(new SPTreeNode()); Cell& new_corner = children[i]->boundary; for(unsigned int d = 0; d < NDims; d++) { - if((i / div) % 2 == 1) + 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]); - div *= 2; } } // Move existing points to correct children - for(unsigned int i = 0; i < size; i++) { - for (auto& child : children) { - if (child->insert(index[i], data, widths, depth+1)) { - break; - } - } + for(unsigned int i = 0; i < QT_NODE_CAPACITY; i++) { + auto c = which_child(data+index[i]*NDims); + children[c]->insert(index[i], data, widths, depth+1); } - - // Empty parent node - size = 0; } template bool SPTreeNode::isCorrect(const double* data, typename vector::const_iterator width) const { - for(unsigned int n = 0; n < size; n++) { - const double* point = data + index[n] * NDims; - if(!boundary.containsPoint(point, *width)) return false; - } - if(!is_leaf()) { + 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)) { @@ -215,11 +223,11 @@ unsigned int SPTreeNode::getAllIndices(unsigned int* indices, unsigned in { // 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()) { + 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); } @@ -243,9 +251,8 @@ void SPTreeNode::computeNonEdgeForces(unsigned int point_index, // Make sure that we spend no time on empty nodes or self-interactions if(cum_size == 0 || - (size == 1 && - __builtin_expect(index[0] == point_index, 0) && - is_leaf())) return; + (cum_size == 1 && + __builtin_expect(index[0] == point_index, 0))) return; // Compute distance between point and center-of-mass double sqdist = .0; @@ -258,8 +265,14 @@ void SPTreeNode::computeNonEdgeForces(unsigned int point_index, // Check whether we can use this node as a "summary" // max_width / sqrt(sqdist) < theta - if(is_leaf() || max_width_squared < theta * theta * sqdist) { - + 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 < theta * theta * sqdist) { // Compute and add t-SNE force between point and current node sqdist = 1.0 / (1.0 + sqdist); double mult = cum_size * sqdist; @@ -268,8 +281,26 @@ void SPTreeNode::computeNonEdgeForces(unsigned int point_index, for(size_t d = 0; d < NDims; ++d) { neg_f[d] += mult * diff[d]; } - } - else { + } else if (QT_NODE_CAPACITY > 1 && cum_size <= QT_NODE_CAPACITY) { + // Need to compute forces on a per-point basis. + auto data = data_point - point_index*NDims; + for (int j = 0; j < cum_size; ++j) { + if (index[j] == point_index) { + continue; + } + auto pt = data+index[j]*NDims; + for(int i = 0; i < NDims; ++i) { + diff[i] = data_point[i] - pt[i]; + sqdist += diff[i] * diff[i]; + } + 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 { // Recursively apply Barnes-Hut to children max_width_squared /= 4.0; @@ -290,11 +321,11 @@ void SPTreeNode::print(const double* data) const if(is_leaf()) { fprintf(stderr,"Leaf node; data = ["); - for(int i = 0; i < size; i++) { + 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 < size - 1) fprintf(stderr,"\n"); + if(i < cum_size - 1) fprintf(stderr,"\n"); else fprintf(stderr,"]\n"); } } @@ -327,7 +358,7 @@ SPTree::SPTree(const double* inp_data, unsigned int N) : data(inp_data) { const double *elem = inp_data; for(unsigned int n = 0; n < N; ++n) { - for(unsigned int d = 0; d < NDims; ++d) { + 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]; @@ -397,10 +428,10 @@ void SPTree::computeEdgeForces(unsigned int* row_P, unsigned int* col_P, // Loop over all edges in the graph unsigned int ind1 = 0; - double sqdist; 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; unsigned int ind2 = col_P[i] * NDims; @@ -431,3 +462,5 @@ void SPTree::print() const { // 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 94000ca..0e5a947 100644 --- a/tsne/bh_sne_src/sptree.h +++ b/tsne/bh_sne_src/sptree.h @@ -38,6 +38,7 @@ #include #include +namespace TSNE { namespace _sptree_internal { template @@ -68,21 +69,19 @@ class SPTreeNode final { 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 - Cell boundary; - // Properties of this node in the tree - unsigned int size; unsigned int cum_size; // Indices in this space-partitioning tree node, corresponding center-of-mass, and list of all children - point_t center_of_mass; std::array index; + // Axis-aligned bounding box stored as a center with half-dimensions to represent the boundaries of this quad tree + Cell boundary; + point_t center_of_mass; + // Children std::array>, no_children> children; @@ -90,6 +89,8 @@ class SPTreeNode final { 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; @@ -155,4 +156,6 @@ class SPTree SPTree(const SPTree&) = delete; }; +} // namespace TSNE + #endif diff --git a/tsne/bh_sne_src/tsne.cpp b/tsne/bh_sne_src/tsne.cpp index 6a3b762..6c1dc2b 100644 --- a/tsne/bh_sne_src/tsne.cpp +++ b/tsne/bh_sne_src/tsne.cpp @@ -92,11 +92,16 @@ void run(double* X, int N, int D, double* Y, double perplexity, double theta, in } // 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); + 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); bool exact = (theta == .0) ? true : false; - 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; @@ -117,6 +122,7 @@ void run(double* X, int N, int D, double* Y, double perplexity, double theta, in for(int i = 0; i < N * D; i++) { if(fabs(X[i]) > max_X) max_X = fabs(X[i]); } + fprintf(stderr,"max deviation from mean == %f\n", max_X); for(int i = 0; i < N * D; i++) X[i] /= max_X; // Compute input similarities for exact t-SNE @@ -125,8 +131,8 @@ void run(double* X, int N, int D, double* Y, double perplexity, double theta, in if(exact) { // Compute similarities - fprintf(stderr,"Exact?"); P.resize(N * N); + fprintf(stderr,"Computing exact perplexity...\n"); computeGaussianPerplexity(X, N, D, P.data(), perplexity); // Symmetrize input similarities @@ -149,6 +155,7 @@ void run(double* X, int N, int D, double* Y, double perplexity, double theta, in // 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)); @@ -642,7 +649,7 @@ void computeSquaredEuclideanDistance(const double* X, int N, int D, double* DD) void zeroMean(double* X, int N, int D) { // Compute data mean - vector mean(0); + vector mean(D, 0); int nD = 0; for(int n = 0; n < N; n++) { for(int d = 0; d < D; d++) { From dedcafe434cdb541fa3b7893815aff64eb38d867 Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Wed, 20 Jun 2018 20:40:50 -0700 Subject: [PATCH 15/22] Trim the fat off VpTree. --- tsne/bh_sne_src/tsne.cpp | 8 ++-- tsne/bh_sne_src/vptree.h | 97 ++++++++++++++++++++-------------------- 2 files changed, 53 insertions(+), 52 deletions(-) diff --git a/tsne/bh_sne_src/tsne.cpp b/tsne/bh_sne_src/tsne.cpp index 6c1dc2b..c8beb89 100644 --- a/tsne/bh_sne_src/tsne.cpp +++ b/tsne/bh_sne_src/tsne.cpp @@ -451,9 +451,9 @@ void computeGaussianPerplexity(const double* X, int N, int D, unsigned int** _ro 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; - vector obj_X(N, DataPoint(D, -1, X)); - for(int n = 0; n < N; n++) obj_X[n] = DataPoint(D, n, X + n * D); + 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 @@ -519,7 +519,7 @@ void computeGaussianPerplexity(const 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(); + col_P[row_P[n] + m] = (unsigned int) (indices[m + 1]._x-X)/D; val_P[row_P[n] + m] = cur_P[m]; } } diff --git a/tsne/bh_sne_src/vptree.h b/tsne/bh_sne_src/vptree.h index a8bd7d7..7e1f819 100644 --- a/tsne/bh_sne_src/vptree.h +++ b/tsne/bh_sne_src/vptree.h @@ -40,52 +40,53 @@ #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 -class DataPoint +namespace TSNE { + +class DataPoint final { - int _ind; public: - std::vector _x; - DataPoint() : _ind(-1) {} - DataPoint(int D, int ind, const double* x) : _ind(ind), _x(x, x+D) {} + const double* _x; + DataPoint() = default; + DataPoint(const double* x) : _x(x) {} DataPoint(const DataPoint& other) = default; - int index() const { return _ind; } - int dimensionality() const { return (int)_x.size(); } double x(int d) const { return _x[d]; } }; -double euclidean_distance(const DataPoint &t1, const DataPoint &t2) { - double dd = .0; - const std::vector& x1 = t1._x; - const std::vector& x2 = t2._x; - double diff; - for(int d = 0; d < x1.size(); 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()); } @@ -101,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(); @@ -117,6 +118,7 @@ class VpTree } private: + const Distance distance; std::vector _items; double _tau; @@ -125,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 @@ -153,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 nullptr; + 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 @@ -181,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]); @@ -219,24 +218,26 @@ class VpTree // 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 From 0ffa88d7f12941e0e22674d4c3b1a04d533367bb Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Thu, 21 Jun 2018 13:06:31 -0700 Subject: [PATCH 16/22] Store child nodes in a contiguous array. Rather than an array of pointers, a pointer to an array. Improved memory locality yeilds a further ~20% performance improvement! --- tsne/bh_sne_src/sptree.cpp | 31 +++++++++++++++++-------------- tsne/bh_sne_src/sptree.h | 2 +- 2 files changed, 18 insertions(+), 15 deletions(-) diff --git a/tsne/bh_sne_src/sptree.cpp b/tsne/bh_sne_src/sptree.cpp index 10547de..00b94cd 100644 --- a/tsne/bh_sne_src/sptree.cpp +++ b/tsne/bh_sne_src/sptree.cpp @@ -154,13 +154,14 @@ bool SPTreeNode::insert(unsigned int new_index, const double* data, vecto cum_size--; return true; } - } // Otherwise, we need to subdivide the current cell + } + // We need to subdivide the current cell subdivide(data, widths, depth); } // Find out where the point can be inserted auto c = which_child(point); - return children[c]->insert(new_index, data, widths, depth+1); + return (*children)[c].insert(new_index, data, widths, depth+1); } // Create four children which fully divide this cell into four quads of equal area @@ -180,9 +181,10 @@ void SPTreeNode::subdivide(const double* data, vector* widths, t const point_t& new_width = (*widths)[depth+1]; // Create new children + children.reset(new array()); + auto& chi = (*children); for(unsigned int i = 0; i < no_children; i++) { - children[i].reset(new SPTreeNode()); - Cell& new_corner = children[i]->boundary; + 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]); @@ -194,7 +196,7 @@ void SPTreeNode::subdivide(const double* data, vector* widths, t // Move existing points to correct children for(unsigned int i = 0; i < QT_NODE_CAPACITY; i++) { auto c = which_child(data+index[i]*NDims); - children[c]->insert(index[i], data, widths, depth+1); + chi[c].insert(index[i], data, widths, depth+1); } } @@ -208,8 +210,8 @@ bool SPTreeNode::isCorrect(const double* data, } } else { ++width; - for(const auto& child : children) { - if (!child->isCorrect(data, width)) { + for(const auto& child : *children) { + if (!child.isCorrect(data, width)) { return false; } } @@ -228,8 +230,8 @@ unsigned int SPTreeNode::getAllIndices(unsigned int* indices, unsigned in loc += cum_size; } else { // Gather indices in children - for(const auto& child : children) - loc = child->getAllIndices(indices, loc); + for(const auto& child : *children) + loc = child.getAllIndices(indices, loc); } return loc; } @@ -238,7 +240,7 @@ 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()); + for(const auto& child : *children) depth = max(depth, child.getDepth()); return 1u + depth; } @@ -281,6 +283,7 @@ void SPTreeNode::computeNonEdgeForces(unsigned int point_index, for(size_t d = 0; d < NDims; ++d) { neg_f[d] += mult * diff[d]; } + /* } else if (QT_NODE_CAPACITY > 1 && cum_size <= QT_NODE_CAPACITY) { // Need to compute forces on a per-point basis. auto data = data_point - point_index*NDims; @@ -299,13 +302,13 @@ void SPTreeNode::computeNonEdgeForces(unsigned int point_index, 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(unsigned int i = 0; i < no_children; i++) - children[i]->computeNonEdgeForces(point_index, data_point, + for(const auto& child : *children) + child.computeNonEdgeForces(point_index, data_point, theta, neg_f, sum_Q, max_width_squared); } } @@ -333,7 +336,7 @@ void SPTreeNode::print(const double* data) const 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(int i = 0; i < no_children; i++) children[i]->print(data); + for(const auto& child : *children) child.print(data); } } diff --git a/tsne/bh_sne_src/sptree.h b/tsne/bh_sne_src/sptree.h index 0e5a947..e6c2836 100644 --- a/tsne/bh_sne_src/sptree.h +++ b/tsne/bh_sne_src/sptree.h @@ -83,7 +83,7 @@ class SPTreeNode final { point_t center_of_mass; // Children - std::array>, no_children> children; + std::unique_ptr, no_children>> children; // Disallow copy SPTreeNode(const SPTreeNode&) = delete; From 8a5efce8d323944463927594f9b0c67bfa4e66e1 Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Thu, 21 Jun 2018 19:32:10 -0700 Subject: [PATCH 17/22] Add alignment directives for better vectorization. Also enable sse3. --- setup.py | 2 +- tsne/bh_sne_src/sptree.cpp | 40 +++++++++++--------------------------- tsne/bh_sne_src/sptree.h | 21 ++++++++++---------- 3 files changed, 23 insertions(+), 40 deletions(-) diff --git a/setup.py b/setup.py index 88b1447..50ad878 100644 --- a/setup.py +++ b/setup.py @@ -45,7 +45,7 @@ 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=['-msse2', '-O3', '-fPIC', '-w', '-ffast-math', '-std=c++11', + 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++')] diff --git a/tsne/bh_sne_src/sptree.cpp b/tsne/bh_sne_src/sptree.cpp index 00b94cd..896c82a 100644 --- a/tsne/bh_sne_src/sptree.cpp +++ b/tsne/bh_sne_src/sptree.cpp @@ -147,8 +147,9 @@ bool SPTreeNode::insert(unsigned int new_index, const double* data, vecto 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] != data[index[n] * NDims + d], 1)) { duplicate = false; break; } + if(__builtin_expect(point[d] != dp[d], 1)) { duplicate = false; break; } } if (__builtin_expect(duplicate, 0)) { cum_size--; @@ -257,9 +258,9 @@ void SPTreeNode::computeNonEdgeForces(unsigned int point_index, __builtin_expect(index[0] == point_index, 0))) return; // Compute distance between point and center-of-mass + alignas(16) point_t diff; double sqdist = .0; - array diff; for(int i = 0; i < NDims; ++i) { diff[i] = data_point[i] - center_of_mass[i]; sqdist += diff[i] * diff[i]; @@ -283,26 +284,6 @@ void SPTreeNode::computeNonEdgeForces(unsigned int point_index, for(size_t d = 0; d < NDims; ++d) { neg_f[d] += mult * diff[d]; } - /* - } else if (QT_NODE_CAPACITY > 1 && cum_size <= QT_NODE_CAPACITY) { - // Need to compute forces on a per-point basis. - auto data = data_point - point_index*NDims; - for (int j = 0; j < cum_size; ++j) { - if (index[j] == point_index) { - continue; - } - auto pt = data+index[j]*NDims; - for(int i = 0; i < NDims; ++i) { - diff[i] = data_point[i] - pt[i]; - sqdist += diff[i] * diff[i]; - } - 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 { // Recursively apply Barnes-Hut to children @@ -354,7 +335,7 @@ 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; - point_t mean_Y, min_Y, max_Y; + 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); @@ -430,18 +411,18 @@ void SPTree::computeEdgeForces(unsigned int* row_P, unsigned int* col_P, { // Loop over all edges in the graph - unsigned int ind1 = 0; + 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; - unsigned int ind2 = col_P[i] * NDims; + const double* data_2 = data + col_P[i] * NDims; - std::array diffs; + alignas(16) point_t diffs; for(unsigned int d = 0; d < NDims; d++) { - diffs[d] = data[ind1 + d] - data[ind2 + d]; + diffs[d] = data_1[d] - data_2[d]; sqdist += diffs[d] * diffs[d]; } @@ -449,10 +430,11 @@ void SPTree::computeEdgeForces(unsigned int* row_P, unsigned int* col_P, // Sum positive force for(unsigned int d = 0; d < NDims; d++) { - pos_f[ind1 + d] += sqdist * diffs[d]; + pos_f[d] += sqdist * diffs[d]; } } - ind1 += NDims; + pos_f += NDims; + data_1 += NDims; } } diff --git a/tsne/bh_sne_src/sptree.h b/tsne/bh_sne_src/sptree.h index e6c2836..990f34a 100644 --- a/tsne/bh_sne_src/sptree.h +++ b/tsne/bh_sne_src/sptree.h @@ -42,7 +42,7 @@ namespace TSNE { namespace _sptree_internal { template -class Cell final { +class alignas(16) Cell final { public: typedef std::array point_t; @@ -62,7 +62,7 @@ class Cell final { }; template -class SPTreeNode final { +class alignas(16) SPTreeNode final { public: typedef typename Cell::point_t point_t; @@ -72,19 +72,19 @@ class SPTreeNode final { // Fixed constants static constexpr unsigned int QT_NODE_CAPACITY = 1; - // 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; - // Axis-aligned bounding box stored as a center with half-dimensions to represent the boundaries of this quad tree - Cell boundary; 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; @@ -115,8 +115,9 @@ class SPTreeNode final { }; template <> -struct SPTreeNode<0> +class SPTreeNode<0> { + public: enum { no_children = 1 }; }; From 2b02650629c070c9add6fbe2fa506365943b3272 Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Thu, 21 Jun 2018 21:26:30 -0700 Subject: [PATCH 18/22] Factor theta into max_width_squared early on. Saves 10% or so. --- tsne/bh_sne_src/sptree.cpp | 10 +++++----- tsne/bh_sne_src/sptree.h | 1 - 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/tsne/bh_sne_src/sptree.cpp b/tsne/bh_sne_src/sptree.cpp index 896c82a..8650f1b 100644 --- a/tsne/bh_sne_src/sptree.cpp +++ b/tsne/bh_sne_src/sptree.cpp @@ -248,7 +248,7 @@ unsigned int SPTreeNode::getDepth() const { // Compute non-edge forces using Barnes-Hut algorithm template void SPTreeNode::computeNonEdgeForces(unsigned int point_index, - const double* data_point, double theta, double neg_f[], double* sum_Q, + const double* data_point, double neg_f[], double* sum_Q, double max_width_squared) const { @@ -275,7 +275,7 @@ void SPTreeNode::computeNonEdgeForces(unsigned int point_index, for(size_t d = 0; d < NDims; ++d) { neg_f[d] += mult * diff[d]; } - } else if (max_width_squared < theta * theta * sqdist) { + } else if (max_width_squared < sqdist) { // Compute and add t-SNE force between point and current node sqdist = 1.0 / (1.0 + sqdist); double mult = cum_size * sqdist; @@ -290,7 +290,7 @@ void SPTreeNode::computeNonEdgeForces(unsigned int point_index, max_width_squared /= 4.0; for(const auto& child : *children) child.computeNonEdgeForces(point_index, data_point, - theta, neg_f, sum_Q, max_width_squared); + neg_f, sum_Q, max_width_squared); } } @@ -401,8 +401,8 @@ unsigned int SPTree::getDepth() const { // 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(); - node.computeNonEdgeForces(point_index, data + point_index * NDims, theta, neg_f, sum_Q, max_width*max_width); + double max_width = maxWidth()/theta; + node.computeNonEdgeForces(point_index, data + point_index * NDims, neg_f, sum_Q, max_width*max_width); } // Computes edge forces diff --git a/tsne/bh_sne_src/sptree.h b/tsne/bh_sne_src/sptree.h index 990f34a..a03ef86 100644 --- a/tsne/bh_sne_src/sptree.h +++ b/tsne/bh_sne_src/sptree.h @@ -105,7 +105,6 @@ class alignas(16) SPTreeNode final { void computeNonEdgeForces(unsigned int point_index, const double* data_point, - double theta, double neg_f[], double* sum_Q, double max_width_squared) const; From aaf77bdbdb4f929573743173f96b8a43c36f5985 Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Mon, 25 Jun 2018 14:16:07 -0700 Subject: [PATCH 19/22] Replace more malloc with vector. Also, tabs -> spaces. --- tsne/bh_sne_src/sptree.cpp | 3 +- tsne/bh_sne_src/tsne.cpp | 302 ++++++++++++++++++------------------- 2 files changed, 151 insertions(+), 154 deletions(-) diff --git a/tsne/bh_sne_src/sptree.cpp b/tsne/bh_sne_src/sptree.cpp index 8650f1b..41959da 100644 --- a/tsne/bh_sne_src/sptree.cpp +++ b/tsne/bh_sne_src/sptree.cpp @@ -267,7 +267,6 @@ void SPTreeNode::computeNonEdgeForces(unsigned int point_index, } // Check whether we can use this node as a "summary" - // max_width / sqrt(sqdist) < theta if(cum_size == 1) { sqdist = 1.0 / (1.0 + sqdist); *sum_Q += sqdist; @@ -276,6 +275,7 @@ void SPTreeNode::computeNonEdgeForces(unsigned int point_index, 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; @@ -285,7 +285,6 @@ void SPTreeNode::computeNonEdgeForces(unsigned int point_index, neg_f[d] += mult * diff[d]; } } else { - // Recursively apply Barnes-Hut to children max_width_squared /= 4.0; for(const auto& child : *children) diff --git a/tsne/bh_sne_src/tsne.cpp b/tsne/bh_sne_src/tsne.cpp index c8beb89..8bba177 100644 --- a/tsne/bh_sne_src/tsne.cpp +++ b/tsne/bh_sne_src/tsne.cpp @@ -41,12 +41,14 @@ #include #include #include +#include #include #include "vptree.h" #include "sptree.h" using std::array; +using std::move; using std::vector; namespace TSNE { @@ -54,7 +56,7 @@ namespace { static inline double sign(double x) { return (x == .0 ? .0 : (x < .0 ? -1.0 : 1.0)); } -void 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); 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 @@ -67,7 +69,9 @@ 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, unsigned int** _row_P, unsigned int** _col_P, double** _val_P, double perplexity, int K); +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); @@ -77,7 +81,7 @@ double randn(); 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 + int max_iter, int stop_lying_iter, int mom_switch_iter ) { // Set random seed @@ -122,12 +126,12 @@ void run(double* X, int N, int D, double* Y, double perplexity, double theta, in for(int i = 0; i < N * D; i++) { if(fabs(X[i]) > max_X) max_X = fabs(X[i]); } - fprintf(stderr,"max deviation from mean == %f\n", max_X); for(int i = 0; i < N * D; i++) X[i] /= max_X; // Compute input similarities for exact t-SNE vector P; - unsigned int* row_P; unsigned int* col_P; double* val_P; + vector row_P, col_P; + vector val_P; if(exact) { // Compute similarities @@ -174,24 +178,24 @@ void run(double* X, int N, int D, double* Y, double perplexity, double theta, in // 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; + Y[i] = randn() * .0001; } } else if (use_init) { for(int i = 0; i < N * no_dims; i++) { - Y[i] = init[i]; + Y[i] = init[i]; } } - // Perform main training loop + // 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.data(), Y, N, dY.data()); - else computeGradient(P.data(), row_P, col_P, val_P, Y, N, dY.data(), theta); + 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); @@ -199,10 +203,10 @@ void run(double* X, int N, int D, double* Y, double perplexity, double theta, in // 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++) Y[i] = Y[i] + uY[i]; // Make solution zero-mean - zeroMean(Y, N); + zeroMean(Y, N); // Stop lying about the P-values after a while, and switch momentum if(iter == stop_lying_iter) { @@ -216,24 +220,18 @@ void run(double* X, int N, int D, double* Y, double perplexity, double theta, in end = clock(); double C = .0; if(exact) C = evaluateError(P.data(), Y, N); - else C = evaluateError(row_P, col_P, val_P, Y, N, theta); // doing approximate computation here! + 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); } - start = clock(); + start = clock(); } } end = clock(); total_time += (float) (end - start) / CLOCKS_PER_SEC; - // Clean up memory - if(!exact) { - 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); } @@ -263,8 +261,8 @@ void computeGradient(double* P, unsigned int* inp_row_P, unsigned int* inp_col_P 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; + // 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 vector DD(N * N); @@ -275,7 +273,7 @@ void computeExactGradient(double* P, const double* Y, int N, double* dC) { double sum_Q = .0; int nN = 0; for(int n = 0; n < N; n++) { - for(int m = 0; m < N; m++) { + for(int m = 0; m < N; m++) { if(n != m) { Q[nN + m] = 1 / (1 + DD[nN + m]); sum_Q += Q[nN + m]; @@ -284,12 +282,12 @@ void computeExactGradient(double* P, const double* Y, int N, 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++) { + 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++) { @@ -297,10 +295,10 @@ void computeExactGradient(double* P, const double* Y, int N, double* dC) { } } mD += D; - } + } nN += N; nD += D; - } + } } @@ -317,7 +315,7 @@ double evaluateError(double* P, const double* Y, int N) { double sum_Q = DBL_MIN; vector Q(N * N); for(int n = 0; n < N; n++) { - for(int m = 0; m < N; m++) { + for(int m = 0; m < N; m++) { if(n != m) { Q[nN + m] = 1 / (1 + DD[nN + m]); sum_Q += Q[nN + m]; @@ -334,7 +332,7 @@ double evaluateError(double* P, const double* Y, int N) { C += P[n] * log((P[n] + FLT_MIN) / (Q[n] + FLT_MIN)); } - return C; + return C; } // Evaluate t-SNE cost function (approximately) @@ -369,83 +367,85 @@ double evaluateError(unsigned int* row_P, unsigned int* col_P, double* val_P, co // Compute input similarities with a fixed perplexity 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 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; + + // 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++; + } + + // Row normalize P + for(int m = 0; m < N; m++) rP[m] /= sum_P; + rP += N; + rD += 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, unsigned int** _row_P, unsigned int** _col_P, double** _val_P, double perplexity, int K) { +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; + _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; @@ -527,12 +527,12 @@ void computeGaussianPerplexity(const double* X, int N, int D, unsigned int** _ro // Symmetrizes a sparse matrix -void 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 vector row_counts(N); @@ -555,10 +555,9 @@ void symmetrizeMatrix(unsigned int** _row_P, unsigned int** _col_P, double** _va 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; @@ -603,9 +602,9 @@ void symmetrizeMatrix(unsigned int** _row_P, unsigned int** _col_P, double** _va 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; + *_row_P = move(sym_row_P); + *_col_P = move(sym_col_P); + *_val_P = move(sym_val_P); } // Compute squared Euclidean distance matrix @@ -648,79 +647,78 @@ void computeSquaredEuclideanDistance(const double* X, int N, int D, double* DD) // Makes data zero-mean void zeroMean(double* X, int N, int D) { - // Compute data mean - vector mean(D, 0); + // 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; - } + } } // Makes data zero-mean template void zeroMean(double* X, int N) { - // Compute data mean - array mean; + // 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]; - } + 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; - } + } } // 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; + 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; } } // namespace // 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 - ) { + 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, From 58c2070891ba91f07c11659d29c485a146bb8da0 Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Mon, 25 Jun 2018 14:26:41 -0700 Subject: [PATCH 20/22] Apply consistent clang-format. --- tsne/bh_sne_src/.clang-format | 3 + tsne/bh_sne_src/sptree.cpp | 279 +++++++++---------- tsne/bh_sne_src/sptree.h | 173 ++++++------ tsne/bh_sne_src/tsne.cpp | 501 ++++++++++++++++++---------------- tsne/bh_sne_src/tsne.h | 12 +- 5 files changed, 510 insertions(+), 458 deletions(-) create mode 100644 tsne/bh_sne_src/.clang-format 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/sptree.cpp b/tsne/bh_sne_src/sptree.cpp index 41959da..4f393c5 100644 --- a/tsne/bh_sne_src/sptree.cpp +++ b/tsne/bh_sne_src/sptree.cpp @@ -46,51 +46,54 @@ using std::vector; namespace TSNE { namespace _sptree_internal { -template +template Cell::Cell(const typename Cell::point_t& p) : corner(p) {} -template +template double Cell::getCorner(unsigned int d) const { return corner[d]; } -template +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; } // Checks whether a point lies in a cell -template -bool Cell::containsPoint(const double* point, const typename Cell::point_t& width) const -{ - 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 (corner[d] - width[d] > point[d]) + return false; + if (corner[d] + width[d] < point[d]) + return false; } return true; } // Constructor for SPTreeNode. -template -SPTreeNode::SPTreeNode(const typename SPTreeNode::point_t& inp_corner) : - boundary(inp_corner), cum_size(0) { +template +SPTreeNode::SPTreeNode( + const typename SPTreeNode::point_t& inp_corner) + : boundary(inp_corner), cum_size(0) { center_of_mass.fill(0.0); } // Constructor for SPTreeNode. -template -SPTreeNode::SPTreeNode() : - cum_size(0) { +template +SPTreeNode::SPTreeNode() : cum_size(0) { center_of_mass.fill(0.0); } // Update the corner position. -template -void SPTreeNode::setCorner(const typename SPTreeNode::point_t& inp_corner) { +template +void SPTreeNode::setCorner( + const typename SPTreeNode::point_t& inp_corner) { boundary.setCorner(inp_corner); } @@ -99,7 +102,7 @@ bool SPTreeNode::is_leaf() const { return NDims == 0 || cum_size <= QT_NODE_CAPACITY; } -template +template unsigned int SPTreeNode::which_child(const double* point) const { unsigned int div = 1; unsigned int i = 0; @@ -112,44 +115,48 @@ unsigned int SPTreeNode::which_child(const double* point) const { return i; } -template -bool SPTreeNode::insert(unsigned int new_index, const double* data, vector* widths, typename vector::size_type depth) -{ +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 const double* point = data + new_index * NDims; - if(!boundary.containsPoint(point, (*widths)[depth])) + if (!boundary.containsPoint(point, (*widths)[depth])) return false; // Online update of cumulative size and center-of-mass unsigned int size = cum_size++; if (size != 0) { - double mult2 = 1.0 / (double) cum_size; - double mult1 = (double) (size) * mult2; + 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]; + 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(size < QT_NODE_CAPACITY) { + if (size < QT_NODE_CAPACITY) { index[size] = new_index; return true; } 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++) { + // 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; } + 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--; @@ -162,32 +169,33 @@ bool SPTreeNode::insert(unsigned int new_index, const double* data, vecto // Find out where the point can be inserted auto c = which_child(point); - return (*children)[c].insert(new_index, data, widths, depth+1); + 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 SPTreeNode::subdivide(const double* data, vector* widths, typename vector::size_type depth) { - +// 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()) { + 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++) { + for (unsigned int d = 0; d < NDims; d++) { child_widths[d] = .5 * width[d]; } } - const point_t& new_width = (*widths)[depth+1]; + const point_t& new_width = (*widths)[depth + 1]; // Create new children children.reset(new array()); auto& chi = (*children); - for(unsigned int i = 0; i < no_children; i++) { + 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) + 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]); @@ -195,23 +203,24 @@ void SPTreeNode::subdivide(const double* data, vector* widths, t } // Move existing points to correct children - 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); + 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); } } -template -bool SPTreeNode::isCorrect(const double* data, - typename vector::const_iterator width) const { +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++) { + for (unsigned int n = 0; n < cum_size; n++) { const double* point = data + index[n] * NDims; - if(!boundary.containsPoint(point, *width)) return false; + if (!boundary.containsPoint(point, *width)) + return false; } } else { ++width; - for(const auto& child : *children) { + for (const auto& child : *children) { if (!child.isCorrect(data, width)) { return false; } @@ -221,57 +230,56 @@ bool SPTreeNode::isCorrect(const double* data, } // Build a list of all indices in SPTree -template -unsigned int SPTreeNode::getAllIndices(unsigned int* indices, unsigned int loc) const -{ - +template +unsigned int SPTreeNode::getAllIndices(unsigned int* indices, + unsigned int loc) const { // Gather indices in current quadrant - if(is_leaf()) { - for(unsigned int i = 0; i < cum_size; i++) indices[loc + i] = index[i]; + 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) + for (const auto& child : *children) loc = child.getAllIndices(indices, loc); } return loc; } -template +template unsigned int SPTreeNode::getDepth() const { - if(is_leaf()) return 1; + if (is_leaf()) + return 1; unsigned int depth = 0; - for(const auto& child : *children) depth = max(depth, child.getDepth()); + for (const auto& child : *children) depth = max(depth, child.getDepth()); return 1u + depth; } // Compute non-edge forces using Barnes-Hut algorithm -template +template void SPTreeNode::computeNonEdgeForces(unsigned int point_index, - const double* data_point, double neg_f[], double* sum_Q, - double max_width_squared) const -{ - + 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 || - (cum_size == 1 && - __builtin_expect(index[0] == point_index, 0))) 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; - for(int i = 0; i < NDims; ++i) { + 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" - if(cum_size == 1) { + 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) { + double mult = sqdist * sqdist; + for (size_t d = 0; d < NDims; ++d) { neg_f[d] += mult * diff[d]; } } else if (max_width_squared < sqdist) { @@ -281,42 +289,42 @@ void SPTreeNode::computeNonEdgeForces(unsigned int point_index, double mult = cum_size * sqdist; *sum_Q += mult; mult *= sqdist; - for(size_t d = 0; d < NDims; ++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); + for (const auto& child : *children) + child.computeNonEdgeForces(point_index, data_point, neg_f, sum_Q, + max_width_squared); } } // Print out tree -template -void SPTreeNode::print(const double* data) const -{ - if(cum_size == 0) { - fprintf(stderr,"Empty node\n"); +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++) { + 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"); + 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); + } 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); } } @@ -324,13 +332,13 @@ void SPTreeNode::print(const double* data) const using namespace _sptree_internal; -template +template double SPTree::maxWidth() const { return *max_element(widths[0].begin(), widths[0].end()); } // Top-node constructor for SPTree -- build tree, too! -template +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; @@ -339,23 +347,25 @@ SPTree::SPTree(const double* inp_data, unsigned int N) : data(inp_data) { 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) { + 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]; + 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; + 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++) { + 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); @@ -363,64 +373,59 @@ SPTree::SPTree(const double* inp_data, unsigned int N) : data(inp_data) { } // Insert a point into the SPTree -template +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); +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 -{ +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 -{ +template +void SPTree::getAllIndices(unsigned int* indices) const { node.getAllIndices(indices, 0); } -template +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); +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 -{ - +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++) { - + 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; const double* data_2 = data + col_P[i] * NDims; alignas(16) point_t diffs; - for(unsigned int d = 0; d < NDims; d++) { + for (unsigned int d = 0; d < NDims; d++) { diffs[d] = data_1[d] - data_2[d]; sqdist += diffs[d] * diffs[d]; } @@ -428,7 +433,7 @@ void SPTree::computeEdgeForces(unsigned int* row_P, unsigned int* col_P, sqdist = val_P[i] / sqdist; // Sum positive force - for(unsigned int d = 0; d < NDims; d++) { + for (unsigned int d = 0; d < NDims; d++) { pos_f[d] += sqdist * diffs[d]; } } @@ -438,7 +443,7 @@ void SPTree::computeEdgeForces(unsigned int* row_P, unsigned int* col_P, } // Print out tree -template +template void SPTree::print() const { node.print(data); } diff --git a/tsne/bh_sne_src/sptree.h b/tsne/bh_sne_src/sptree.h index a03ef86..b2ee840 100644 --- a/tsne/bh_sne_src/sptree.h +++ b/tsne/bh_sne_src/sptree.h @@ -41,119 +41,122 @@ namespace TSNE { namespace _sptree_internal { -template +template class alignas(16) Cell final { public: - typedef std::array point_t; + typedef std::array point_t; - Cell() = default; - explicit Cell(const point_t&); + 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; + 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; private: - point_t corner; + point_t corner; - // disallow copy - Cell(const Cell&) = delete; + // disallow copy + Cell(const Cell&) = delete; }; template class 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; + public: + typedef typename Cell::point_t point_t; + enum { no_children = 2 * SPTreeNode::no_children }; - // Properties of this node in the tree - unsigned int cum_size; + private: + // Fixed constants + static constexpr unsigned int QT_NODE_CAPACITY = 1; - // Indices in this space-partitioning tree node, corresponding center-of-mass, and list of all children - std::array index; + // 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; - // Disallow copy - SPTreeNode(const SPTreeNode&) = delete; + // Children + std::unique_ptr, no_children>> children; - 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); + // Properties of this node in the tree + unsigned int cum_size; - bool is_leaf() const; + // Indices in this space-partitioning tree node, corresponding center-of-mass, + // and list of all children + std::array index; -public: - SPTreeNode(); - explicit SPTreeNode(const point_t& corner); + // Disallow copy + SPTreeNode(const SPTreeNode&) = delete; - void setCorner(const point_t& corner); + 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 insert(unsigned int new_index, const double* data, std::vector* widths, typename std::vector::size_type depth); - bool isCorrect(const double* data, typename std::vector::const_iterator width) const; + bool is_leaf() const; - void computeNonEdgeForces(unsigned int point_index, - const double* data_point, - double neg_f[], - double* sum_Q, - double max_width_squared) const; - unsigned int getAllIndices(unsigned int* indices, unsigned int loc) const; - unsigned int getDepth() const; - void print(const double* data) 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); + bool isCorrect(const double* data, + typename std::vector::const_iterator width) const; + + void computeNonEdgeForces(unsigned int point_index, const double* data_point, + double neg_f[], double* sum_Q, + double max_width_squared) const; + unsigned int getAllIndices(unsigned int* indices, unsigned int loc) const; + unsigned int getDepth() const; + void print(const double* data) const; }; template <> -class SPTreeNode<0> -{ +class SPTreeNode<0> { public: - enum { no_children = 1 }; + enum { no_children = 1 }; }; } // namespace _sptree_internal -template -class 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); - - bool isCorrect() const; - void getAllIndices(unsigned int* indices) const; - unsigned int getDepth() const; - void computeNonEdgeForces(unsigned int point_index, double theta, double neg_f[], double* sum_Q) const; - void computeEdgeForces(unsigned int* row_P, unsigned int* col_P, double* val_P, int N, double* pos_f) const; - void print() const; - -private: - void fill(unsigned int N); - // Disallow copy - SPTree(const SPTree&) = delete; +template +class 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); + + bool isCorrect() const; + void getAllIndices(unsigned int* indices) const; + unsigned int getDepth() const; + void computeNonEdgeForces(unsigned int point_index, double theta, + double neg_f[], double* sum_Q) const; + void computeEdgeForces(unsigned int* row_P, unsigned int* col_P, + double* val_P, int N, double* pos_f) const; + void print() const; + + private: + void fill(unsigned int N); + // Disallow copy + SPTree(const SPTree&) = delete; }; } // namespace TSNE diff --git a/tsne/bh_sne_src/tsne.cpp b/tsne/bh_sne_src/tsne.cpp index 8bba177..12dbadf 100644 --- a/tsne/bh_sne_src/tsne.cpp +++ b/tsne/bh_sne_src/tsne.cpp @@ -44,8 +44,8 @@ #include #include -#include "vptree.h" #include "sptree.h" +#include "vptree.h" using std::array; using std::move; @@ -54,24 +54,33 @@ using std::vector; namespace TSNE { namespace { -static inline double sign(double x) { return (x == .0 ? .0 : (x < .0 ? -1.0 : 1.0)); } +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); +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); +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); +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, 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); + 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); @@ -79,32 +88,33 @@ 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 - ) { - +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"); + 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); + fprintf(stderr, + "Using D = %d, no_dims = %d, perplexity = %f, and theta = %f\n", D, + no_dims, perplexity, theta); bool exact = (theta == .0) ? true : false; - fprintf(stderr,"Using max_iter = %d, stop_lying_iter = %d, mom_switch_iter = %d\n", + 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 @@ -119,127 +129,158 @@ void run(double* X, int N, int D, double* Y, double perplexity, double theta, in 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 vector P; vector row_P, col_P; vector val_P; - if(exact) { - + if (exact) { // Compute similarities P.resize(N * N); - fprintf(stderr,"Computing exact perplexity...\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"); + 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)); + 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.data(), Y, N, dY.data()); - else computeGradient(P.data(), row_P.data(), col_P.data(), val_P.data(), Y, N, dY.data(), 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); // 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.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); + 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(); } } - end = clock(); total_time += (float) (end - start) / CLOCKS_PER_SEC; + end = clock(); + total_time += (float)(end - start) / CLOCKS_PER_SEC; - 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) 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) -{ +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(Y, N); @@ -249,10 +290,11 @@ void computeGradient(double* P, unsigned int* inp_row_P, unsigned int* inp_col_P 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); + 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 < len; i++) { + for (int i = 0; i < len; i++) { dC[i] = pos_f[i] - (neg_f[i] / sum_Q); } } @@ -260,9 +302,8 @@ void computeGradient(double* P, unsigned int* inp_row_P, unsigned int* inp_col_P // Compute gradient of the t-SNE cost function (exact) 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; + for (int i = 0; i < N * D; i++) dC[i] = 0.0; // Compute the squared Euclidean distance matrix vector DD(N * N); @@ -272,9 +313,9 @@ void computeExactGradient(double* P, const double* Y, int N, double* dC) { 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]; } @@ -285,12 +326,12 @@ void computeExactGradient(double* P, const double* Y, int N, double* dC) { // 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; } } @@ -301,11 +342,9 @@ void computeExactGradient(double* P, const double* Y, int N, double* dC) { } } - // Evaluate t-SNE cost function (exactly) template double evaluateError(double* P, const double* Y, int N) { - // Compute the squared Euclidean distance matrix vector DD(N * N); computeSquaredEuclideanDistance(Y, N, DD.data()); @@ -314,22 +353,22 @@ double evaluateError(double* P, const double* Y, int N) { int nN = 0; double sum_Q = DBL_MIN; vector Q(N * N); - 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]; - } - 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)); } return C; @@ -337,25 +376,26 @@ double evaluateError(double* P, const double* Y, int N) { // Evaluate t-SNE cost function (approximately) template -double evaluateError(unsigned int* row_P, unsigned int* col_P, double* val_P, const double* Y, int N, double theta) -{ +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(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)); } @@ -363,59 +403,54 @@ double evaluateError(unsigned int* row_P, unsigned int* col_P, double* val_P, co return C; } - // Compute input similarities with a fixed perplexity -void computeGaussianPerplexity(const double* X, int N, int D, double* P, double perplexity) { - +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 - double *rP = P; - double *rD = DD.data(); - for(int n = 0; n < N; n++) { - + 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 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) { - + while (!found && iter < 200) { // Compute Gaussian kernel row - for(int m = 0; m < N; m++) rP[m] = exp(-beta * rD[m]); + for (int m = 0; m < N; m++) rP[m] = exp(-beta * rD[m]); rP[n] = DBL_MIN; // Compute entropy of current row sum_P = DBL_MIN; - for(int m = 0; m < N; m++) sum_P += rP[m]; + 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]; + 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) { + 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; @@ -427,42 +462,46 @@ void computeGaussianPerplexity(const double* X, int N, int D, double* P, double } // Row normalize P - for(int m = 0; m < N; m++) rP[m] /= sum_P; + for (int m = 0; m < N; m++) rP[m] /= sum_P; rP += N; rD += 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"); +// 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->resize(N+1); - _col_P->resize(N*K,0); - _val_P->resize(N*K,0); + _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((euclidean_distance(D))); + VpTree tree((euclidean_distance(D))); vector obj_X(N); - for(int n = 0; n < N; n++) obj_X[n] = DataPoint(X + n * D); + 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(); @@ -473,39 +512,39 @@ void computeGaussianPerplexity(const double* X, int N, int D, vector 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; @@ -517,18 +556,18 @@ void computeGaussianPerplexity(const double* X, int N, int D, vector* _row_P, vector* _col_P, vector* _val_P, int N) { - +void symmetrizeMatrix(vector* _row_P, + vector* _col_P, vector* _val_P, + int N) { // Get sparse matrix vector& row_P = *_row_P; vector& col_P = *_col_P; @@ -536,15 +575,16 @@ void symmetrizeMatrix(vector* _row_P, vector* _col_P // Count number of elements and row counts of symmetric matrix vector row_counts(N); - for(int n = 0; n < N; n++) { - for(int i = row_P[n]; i < row_P[n + 1]; i++) { - + 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]]++; @@ -552,7 +592,7 @@ void symmetrizeMatrix(vector* _row_P, vector* _col_P } } 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 vector sym_row_P(N + 1); @@ -561,45 +601,52 @@ void symmetrizeMatrix(vector* _row_P, vector* _col_P // 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 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]) + 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 *_row_P = move(sym_row_P); @@ -611,14 +658,14 @@ void symmetrizeMatrix(vector* _row_P, vector* _col_P 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; @@ -626,16 +673,17 @@ void computeSquaredEuclideanDistance(const double* X, int N, double* DD) { } } -void computeSquaredEuclideanDistance(const 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) { + 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; @@ -643,27 +691,25 @@ void computeSquaredEuclideanDistance(const double* X, int N, int D, double* DD) } } - // Makes data zero-mean 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++) { + 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 nD = 0; - for(int n = 0; n < N; n++) { - for(int d = 0; d < D; d++) { + for (int n = 0; n < N; n++) { + for (int d = 0; d < D; d++) { X[nD + d] -= mean[d]; } nD += D; @@ -673,40 +719,38 @@ void zeroMean(double* X, int N, int D) { // 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++) { + 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 nD = 0; - for(int n = 0; n < N; n++) { - for(int d = 0; d < D; d++) { + for (int n = 0; n < N; n++) { + for (int d = 0; d < D; d++) { X[nD + d] -= mean[d]; } nD += D; } } - // 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; + 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)); + } while ((radius >= 1.0) || (radius == 0.0)); radius = sqrt(-2 * log(radius) / radius); x *= radius; y *= radius; @@ -716,22 +760,21 @@ double randn() { } // namespace // 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"); +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"); } } diff --git a/tsne/bh_sne_src/tsne.h b/tsne/bh_sne_src/tsne.h index caa4fbf..27aca03 100644 --- a/tsne/bh_sne_src/tsne.h +++ b/tsne/bh_sne_src/tsne.h @@ -30,20 +30,18 @@ * */ - #ifndef TSNE_H #define TSNE_H namespace 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=1000, - int stop_lying_iter=250, int mom_switch_iter=250 - ); + 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); + 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); From c6443027df4f7f73554fe56ff5b3b542f425307d Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Wed, 27 Jun 2018 09:58:49 -0700 Subject: [PATCH 21/22] Factor beta multiplication. Why do it N times when you can do it just once? --- tsne/bh_sne_src/tsne.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tsne/bh_sne_src/tsne.cpp b/tsne/bh_sne_src/tsne.cpp index 12dbadf..ba9ed9b 100644 --- a/tsne/bh_sne_src/tsne.cpp +++ b/tsne/bh_sne_src/tsne.cpp @@ -528,8 +528,8 @@ void computeGaussianPerplexity(const double* X, int N, int D, 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); + 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); From 90c3647f138bc227ccbffd354539f50ca981ab53 Mon Sep 17 00:00:00 2001 From: Adam Azarchs Date: Wed, 27 Jun 2018 14:07:47 -0700 Subject: [PATCH 22/22] Add optimizer hints for hot paths. Also add visibility attributes for internal classes. This allows for more aggressive linker garbage collection and optimization. Turning on -fvisiblity=hidden globablly would require cython to set visibility attributes on the entry points it exports, so we're not going to get the full benefit, unfortunately. Also change containsPoint comparison to abs(corner-point)>width. This replaces a second subtract/compare with a single vectorized bitwise and instruction, and by reducing branch counts also helps CPU pipelining. --- tsne/bh_sne_src/sptree.cpp | 6 ++--- tsne/bh_sne_src/sptree.h | 51 +++++++++++++++++++++----------------- tsne/bh_sne_src/tsne.cpp | 3 +-- 3 files changed, 32 insertions(+), 28 deletions(-) diff --git a/tsne/bh_sne_src/sptree.cpp b/tsne/bh_sne_src/sptree.cpp index 4f393c5..40f6f1a 100644 --- a/tsne/bh_sne_src/sptree.cpp +++ b/tsne/bh_sne_src/sptree.cpp @@ -34,9 +34,11 @@ #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; @@ -69,9 +71,7 @@ template bool Cell::containsPoint( const double* point, const typename Cell::point_t& width) const { for (int d = 0; d < NDims; ++d) { - if (corner[d] - width[d] > point[d]) - return false; - if (corner[d] + width[d] < point[d]) + if (abs(corner[d] - point[d]) > width[d]) return false; } return true; diff --git a/tsne/bh_sne_src/sptree.h b/tsne/bh_sne_src/sptree.h index b2ee840..ba7b325 100644 --- a/tsne/bh_sne_src/sptree.h +++ b/tsne/bh_sne_src/sptree.h @@ -30,7 +30,6 @@ * */ - #ifndef SPTREE_H #define SPTREE_H @@ -42,7 +41,7 @@ namespace TSNE { namespace _sptree_internal { template -class alignas(16) Cell final { +class[[gnu::visibility("internal")]] alignas(16) Cell final { public: typedef std::array point_t; @@ -62,7 +61,7 @@ class alignas(16) Cell final { }; template -class alignas(16) SPTreeNode final { +class[[gnu::visibility("internal")]] alignas(16) SPTreeNode final { public: typedef typename Cell::point_t point_t; enum { no_children = 2 * SPTreeNode::no_children }; @@ -82,8 +81,8 @@ class alignas(16) SPTreeNode final { // 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 + // Indices in this space-partitioning tree node, corresponding + // center-of-mass, and list of all children std::array index; // Disallow copy @@ -105,15 +104,19 @@ class alignas(16) SPTreeNode final { bool insert(unsigned int new_index, const double* data, std::vector* widths, typename std::vector::size_type depth); - bool isCorrect(const double* data, - typename std::vector::const_iterator width) const; - - void computeNonEdgeForces(unsigned int point_index, const double* data_point, - double neg_f[], double* sum_Q, - double max_width_squared) const; - unsigned int getAllIndices(unsigned int* indices, unsigned int loc) const; - unsigned int getDepth() const; - void print(const double* data) const; + 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 <> @@ -125,7 +128,7 @@ class SPTreeNode<0> { } // namespace _sptree_internal template -class SPTree { +class [[gnu::visibility("internal")]] SPTree { public: typedef typename _sptree_internal::SPTreeNode::point_t point_t; enum { no_children = _sptree_internal::SPTreeNode::no_children }; @@ -144,14 +147,16 @@ class SPTree { public: SPTree(const double* inp_data, unsigned int N); - bool isCorrect() const; - void getAllIndices(unsigned int* indices) const; - unsigned int getDepth() const; - void computeNonEdgeForces(unsigned int point_index, double theta, - double neg_f[], double* sum_Q) const; - void computeEdgeForces(unsigned int* row_P, unsigned int* col_P, - double* val_P, int N, double* pos_f) const; - void print() const; + 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); diff --git a/tsne/bh_sne_src/tsne.cpp b/tsne/bh_sne_src/tsne.cpp index ba9ed9b..8f2368e 100644 --- a/tsne/bh_sne_src/tsne.cpp +++ b/tsne/bh_sne_src/tsne.cpp @@ -30,7 +30,6 @@ * */ - #include "tsne.h" #include @@ -111,7 +110,7 @@ void run(double* X, int N, int D, double* Y, double perplexity, double theta, fprintf(stderr, "Using D = %d, no_dims = %d, perplexity = %f, and theta = %f\n", D, no_dims, perplexity, theta); - bool exact = (theta == .0) ? true : false; + 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",