Rename ndim_ to num_dim_ for more clarity and readability.

Also renamed the initialize() overload arguments for clarity.

include/qdtsne/SPTree.hpp

@@ -42,7 +42,7 @@ namespace qdtsne {
 
 namespace internal {
 
-template<int ndim_, typename Float_>
+template<int num_dim_, typename Float_>
 class SPTree {
 public:
     SPTree(size_t npts, int maxdepth) : my_npts(npts), my_maxdepth(maxdepth), my_locations(my_npts) {
@@ -53,7 +53,7 @@ class SPTree {
 public:
     struct Node {
         Node(size_t i, const Float_* point) : index(i) {
-            std::copy_n(point, ndim_, center_of_mass.data());
+            std::copy_n(point, num_dim_, center_of_mass.data());
             fill();
             return;
         }
@@ -65,11 +65,11 @@ class SPTree {
         }
 
     public:
-        static constexpr int nchildren = (1 << ndim_); 
+        static constexpr int nchildren = (1 << num_dim_); 
 
         std::array<size_t, nchildren> children;
-        std::array<Float_, ndim_> midpoint, halfwidth;
-        std::array<Float_, ndim_> center_of_mass;
+        std::array<Float_, num_dim_> midpoint, halfwidth;
+        std::array<Float_, num_dim_> center_of_mass;
         Float_ max_width = 0;
 
         size_t number = 1;
@@ -83,8 +83,8 @@ class SPTree {
 
     private:
         void fill() {
-            std::fill_n(midpoint.begin(), ndim_, 0);
-            std::fill_n(halfwidth.begin(), ndim_, 0);
+            std::fill_n(midpoint.begin(), num_dim_, 0);
+            std::fill_n(halfwidth.begin(), num_dim_, 0);
             std::fill_n(children.begin(), nchildren, 0);
         }
     };
@@ -115,39 +115,39 @@ class SPTree {
             my_store[0].is_leaf = false;
             my_store[0].number = my_npts;
 
-            std::array<Float_, ndim_> min_Y{}, max_Y{};
-            std::fill_n(min_Y.begin(), ndim_, std::numeric_limits<Float_>::max());
-            std::fill_n(max_Y.begin(), ndim_, std::numeric_limits<Float_>::lowest());
+            std::array<Float_, num_dim_> min_Y{}, max_Y{};
+            std::fill_n(min_Y.begin(), num_dim_, std::numeric_limits<Float_>::max());
+            std::fill_n(max_Y.begin(), num_dim_, std::numeric_limits<Float_>::lowest());
 
             // Setting the initial midpoint to the center of mass of all points
             // so that the partitioning effectively divides up the points.
             auto& mean_Y = my_store[0].midpoint;
             auto copy = Y;
             for (size_t n = 0; n < my_npts; ++n) {
-                for (int d = 0; d < ndim_; ++d) {
+                for (int d = 0; d < num_dim_; ++d) {
                     auto curval = copy[d];
                     mean_Y[d] += curval;
                     min_Y[d] = std::min(min_Y[d], curval);
                     max_Y[d] = std::max(max_Y[d], curval);
                 }
-                copy += ndim_;
+                copy += num_dim_;
             }
 
-            for (int d = 0; d < ndim_; ++d) {
+            for (int d = 0; d < num_dim_; ++d) {
                 mean_Y[d] /= my_npts;
             }
 
             auto& halfwidth = my_store[0].halfwidth;
-            for (int d = 0; d < ndim_; ++d) {
+            for (int d = 0; d < num_dim_; ++d) {
                 auto mean = mean_Y[d];
                 halfwidth[d] = std::max(max_Y[d] - mean, mean - min_Y[d]) + static_cast<Float_>(1e-5);
             }
         }
 
         auto point = Y;
         my_first_assignment.resize(my_npts);
-        for (size_t i = 0; i < my_npts; ++i, point += ndim_) {
-            std::array<bool, ndim_> side;
+        for (size_t i = 0; i < my_npts; ++i, point += num_dim_) {
+            std::array<bool, num_dim_> side;
             size_t parent = 0;
 
             for (int depth = 1; depth <= my_maxdepth; ++depth) {
@@ -172,10 +172,10 @@ class SPTree {
                     // No need to update the center of mass because it's the same point.
                     int nsame = 0;
                     const auto& center = my_store[current_loc].center_of_mass;
-                    for (int d = 0; d < ndim_; ++d) {
+                    for (int d = 0; d < num_dim_; ++d) {
                         nsame += (center[d] == point[d]);
                     }
-                    if (nsame == ndim_) {
+                    if (nsame == num_dim_) {
                         ++(my_store[current_loc].number);
                         my_first_assignment[i] = my_store[current_loc].index;
                         break;
@@ -209,7 +209,7 @@ class SPTree {
                 const Float_ cum_size = node.number;
                 const Float_ mult1 = (cum_size - 1) / cum_size;
 
-                for (int d = 0; d < ndim_; ++d) {
+                for (int d = 0; d < num_dim_; ++d) {
                     node.center_of_mass[d] *= mult1;
                     node.center_of_mass[d] += point[d] / cum_size;
                 }
@@ -240,7 +240,7 @@ class SPTree {
     size_t find_child (size_t parent, const Float_* point, bool * side) const {
         int multiplier = 1;
         size_t child = 0;
-        for (int d = 0; d < ndim_; ++d) {
+        for (int d = 0; d < num_dim_; ++d) {
             side[d] = (point[d] >= my_store[parent].midpoint[d]);
             child += multiplier * side[d];
             multiplier *= 2;
@@ -251,7 +251,7 @@ class SPTree {
     void set_child_boundaries(size_t parent, size_t child, const bool* keep) {
         auto& current = my_store[child];
         auto& parental = my_store[parent];
-        for (int d = 0; d < ndim_; ++d) {
+        for (int d = 0; d < num_dim_; ++d) {
             current.halfwidth[d] = parental.halfwidth[d] / static_cast<Float_>(2);
             if (keep[d]) {
                 current.midpoint[d] = parental.midpoint[d] + current.halfwidth[d];
@@ -272,43 +272,43 @@ class SPTree {
      *** Non-edge force calculations ***
      ***********************************/
 private:
-    static Float_ compute_sqdist(const Float_* point, const std::array<Float_, ndim_>& center) {
+    static Float_ compute_sqdist(const Float_* point, const std::array<Float_, num_dim_>& center) {
         Float_ sqdist = 0;
-        for (int d = 0; d < ndim_; ++d) {
+        for (int d = 0; d < num_dim_; ++d) {
             Float_ delta = point[d] - center[d];
             sqdist += delta * delta;
         }
         return sqdist;
     }
 
-    static void add_non_edge_forces(const Float_* point, const std::array<Float_, ndim_>& center, Float_ sqdist, size_t count, Float_& result_sum, Float_* neg_f) {
+    static void add_non_edge_forces(const Float_* point, const std::array<Float_, num_dim_>& center, Float_ sqdist, size_t count, Float_& result_sum, Float_* neg_f) {
         const Float_ div = static_cast<Float_>(1) / (static_cast<Float_>(1) + sqdist);
         Float_ mult = count * div;
         result_sum += mult;
         mult *= div;
 #ifdef _OPENMP
         #pragma omp simd
 #endif
-        for (int d = 0; d < ndim_; ++d) {
+        for (int d = 0; d < num_dim_; ++d) {
             neg_f[d] += mult * (point[d] - center[d]);
         }
     }
 
-    static void remove_self_from_center(const Float_* point, const std::array<Float_, ndim_>& center, Float_ count, std::array<Float_, ndim_>& temp) {
+    static void remove_self_from_center(const Float_* point, const std::array<Float_, num_dim_>& center, Float_ count, std::array<Float_, num_dim_>& temp) {
 #ifdef _OPENMP
         #pragma omp simd
 #endif
-        for (int d = 0; d < ndim_; ++d) { 
+        for (int d = 0; d < num_dim_; ++d) { 
             temp[d] = (center[d] * count - point[d]) / (count - 1);
         }
     }
 
 public:
     Float_ compute_non_edge_forces(size_t index, Float_ theta, Float_* neg_f) const {
         Float_ result_sum = 0;
-        const Float_ * point = my_data + index * static_cast<size_t>(ndim_); // cast to avoid overflow.
+        const Float_ * point = my_data + index * static_cast<size_t>(num_dim_); // cast to avoid overflow.
         const auto& cur_children = my_store[0].children;
-        std::fill_n(neg_f, ndim_, 0);
+        std::fill_n(neg_f, num_dim_, 0);
 
         for (int i = 0; i < Node::nchildren; ++i) {
             if (cur_children[i]) {
@@ -323,7 +323,7 @@ class SPTree {
     Float_ compute_non_edge_forces(size_t index, const Float_* point, Float_ theta, Float_* neg_f, size_t position) const {
         const auto& node = my_store[position];
 
-        std::array<Float_, ndim_> temp;
+        std::array<Float_, num_dim_> temp;
         auto center = &(node.center_of_mass);
         size_t count = node.number;
 
@@ -366,7 +366,7 @@ class SPTree {
 public:
     struct LeafApproxWorkspace {
         std::vector<size_t> leaf_indices;
-        std::vector<std::array<Float_, ndim_> > leaf_neg_f;
+        std::vector<std::array<Float_, num_dim_> > leaf_neg_f;
         std::vector<Float_> leaf_sums;
     };
 
@@ -378,7 +378,7 @@ class SPTree {
         auto process_leaf_node = [&](size_t leaf) {
             Float_ result_sum = 0;
             auto neg_f = workspace.leaf_neg_f[leaf].data();
-            std::fill_n(neg_f, ndim_, 0);
+            std::fill_n(neg_f, num_dim_, 0);
 
             const auto& cur_children = my_store[0].children;
             for (int i = 0; i < Node::nchildren; ++i) {
@@ -445,8 +445,8 @@ class SPTree {
 
         const auto& node = my_store[node_loc];
         if (node.number != 1) {
-            const Float_ * point = my_data + index * static_cast<size_t>(ndim_); // cast to avoid overflow.
-            std::array<Float_, ndim_> temp;
+            const Float_ * point = my_data + index * static_cast<size_t>(num_dim_); // cast to avoid overflow.
+            std::array<Float_, num_dim_> temp;
             remove_self_from_center(point, node.center_of_mass, node.number, temp);
             Float_ sqdist = compute_sqdist(point, temp);
             add_non_edge_forces(point, temp, sqdist, node.number - 1, result_sum, neg_f);

include/qdtsne/Status.hpp

@@ -53,26 +53,26 @@ namespace qdtsne {
 /**
  * @brief Current status of the t-SNE iterations.
  *
- * @tparam ndim_ Number of dimensions in the t-SNE embedding.
+ * @tparam num_dim_ Number of dimensions in the t-SNE embedding.
  * @tparam Index_ Integer type for the neighbor indices.
  * @tparam Float_ Floating-point type for the distances.
  *
  * This class holds the precomputed structures required to perform the t-SNE iterations.
  * Instances should not be constructed directly but instead created by `initialize()`.
  */
-template<int ndim_, typename Index_, typename Float_> 
+template<int num_dim_, typename Index_, typename Float_> 
 class Status {
 public:
     /**
      * @cond
      */
     Status(NeighborList<Index_, Float_> neighbors, Options options) :
         my_neighbors(std::move(neighbors)),
-        my_dY(my_neighbors.size() * ndim_), 
-        my_uY(my_neighbors.size() * ndim_), 
-        my_gains(my_neighbors.size() * ndim_, 1.0), 
-        my_pos_f(my_neighbors.size() * ndim_), 
-        my_neg_f(my_neighbors.size() * ndim_), 
+        my_dY(my_neighbors.size() * num_dim_), 
+        my_uY(my_neighbors.size() * num_dim_), 
+        my_gains(my_neighbors.size() * num_dim_, 1.0), 
+        my_pos_f(my_neighbors.size() * num_dim_), 
+        my_neg_f(my_neighbors.size() * num_dim_), 
         my_tree(my_neighbors.size(), options.max_depth),
         my_options(std::move(options))
     {
@@ -88,7 +88,7 @@ class Status {
     NeighborList<Index_, Float_> my_neighbors; 
     std::vector<Float_> my_dY, my_uY, my_gains, my_pos_f, my_neg_f;
 
-    internal::SPTree<ndim_, Float_> my_tree;
+    internal::SPTree<num_dim_, Float_> my_tree;
     std::vector<Float_> my_parallel_buffer; // Buffer to hold parallel-computed results prior to reduction.
 
     Options my_options;
@@ -136,7 +136,7 @@ class Status {
      * Run the algorithm to the specified number of iterations.
      * This can be invoked repeatedly with increasing `limit` to run the algorithm incrementally.
      *
-     * @param[in, out] Y Pointer to a array containing a column-major matrix with number of rows and columns equal to `ndim_` and `num_observations()`, respectively.
+     * @param[in, out] Y Pointer to a array containing a column-major matrix with number of rows and columns equal to `num_dim_` and `num_observations()`, respectively.
      * Each row corresponds to a dimension of the embedding while each column corresponds to an observation.
      * On input, this should contain the initial location of each observation; on output, it is updated to the t-SNE location at the specified number of iterations.
      * @param limit Number of iterations to run up to.
@@ -165,7 +165,7 @@ class Status {
      * If `run()` has already been invoked with an iteration limit, this method will only perform the remaining iterations required for `iteration()` to reach `max_iter()`.
      * If `iteration()` is already greater than `max_iter()`, this method is a no-op.
      *
-     * @param[in, out] Y Pointer to a array containing a column-major matrix with number of rows and columns equal to `ndim_` and `num_observations()`, respectively.
+     * @param[in, out] Y Pointer to a array containing a column-major matrix with number of rows and columns equal to `num_dim_` and `num_observations()`, respectively.
      * Each row corresponds to a dimension of the embedding while each column corresponds to an observation.
      * On input, this should contain the initial location of each observation; on output, it is updated to the t-SNE location at the specified number of iterations.
      */
@@ -207,18 +207,18 @@ class Status {
 
         // Make solution zero-mean
         size_t N = num_observations();
-        for (int d = 0; d < ndim_; ++d) {
+        for (int d = 0; d < num_dim_; ++d) {
             auto start = Y + d;
 
             // Compute means from column-major coordinates.
             Float_ sum = 0;
-            for (size_t i = 0; i < N; ++i, start += ndim_) {
+            for (size_t i = 0; i < N; ++i, start += num_dim_) {
                 sum += *start;
             }
             sum /= N;
 
             start = Y + d;
-            for (size_t i = 0; i < N; ++i, start += ndim_) {
+            for (size_t i = 0; i < N; ++i, start += num_dim_) {
                 *start -= sum;
             }
         }
@@ -237,7 +237,7 @@ class Status {
         Float_ sum_Q = compute_non_edge_forces();
 
         // Compute final t-SNE gradient
-        size_t ntotal = N * static_cast<size_t>(ndim_);
+        size_t ntotal = N * static_cast<size_t>(num_dim_);
         for (size_t i = 0; i < ntotal; ++i) {
             my_dY[i] = my_pos_f[i] - (my_neg_f[i] / sum_Q);
         }
@@ -262,20 +262,20 @@ class Status {
 #endif
 
                 const auto& current = my_neighbors[n];
-                size_t offset = n * static_cast<size_t>(ndim_); // cast to avoid overflow.
+                size_t offset = n * static_cast<size_t>(num_dim_); // cast to avoid overflow.
                 const Float_* self = Y + offset;
                 Float_* pos_out = my_pos_f.data() + offset;
 
                 for (const auto& x : current) {
                     Float_ sqdist = 0; 
-                    const Float_* neighbor = Y + static_cast<size_t>(x.first) * ndim_; // cast to avoid overflow.
-                    for (int d = 0; d < ndim_; ++d) {
+                    const Float_* neighbor = Y + static_cast<size_t>(x.first) * num_dim_; // cast to avoid overflow.
+                    for (int d = 0; d < num_dim_; ++d) {
                         Float_ delta = self[d] - neighbor[d];
                         sqdist += delta * delta;
                     }
 
                     const Float_ mult = multiplier * x.second / (static_cast<Float_>(1) + sqdist);
-                    for (int d = 0; d < ndim_; ++d) {
+                    for (int d = 0; d < num_dim_; ++d) {
                         pos_out[d] += mult * (self[d] - neighbor[d]);
                     }
                 }
@@ -317,7 +317,7 @@ class Status {
                 for (size_t n = first_; n < last_; ++n) {
 #endif                
 
-                    auto neg_ptr = my_neg_f.data() + n * static_cast<size_t>(ndim_); // cast to avoid overflow.
+                    auto neg_ptr = my_neg_f.data() + n * static_cast<size_t>(num_dim_); // cast to avoid overflow.
                     if (my_options.leaf_approximation) {
                         my_parallel_buffer[n] = my_tree.compute_non_edge_forces_from_leaves(n, neg_ptr, my_leaf_workspace);
                     } else {
@@ -338,7 +338,7 @@ class Status {
 
         Float_ sum_Q = 0;
         for (size_t n = 0; n < N; ++n) {
-            auto neg_ptr = my_neg_f.data() + n * static_cast<size_t>(ndim_); // cast to avoid overflow.
+            auto neg_ptr = my_neg_f.data() + n * static_cast<size_t>(num_dim_); // cast to avoid overflow.
             if (my_options.leaf_approximation) {
                 sum_Q += my_tree.compute_non_edge_forces_from_leaves(n, neg_ptr, my_leaf_workspace);
             } else {