complete!

Author: jinzx10

source/module_base/grid/batch.cpp

@@ -24,22 +24,21 @@ namespace {
  *
  * It can be shown that the optimal c is the centroid of the points, and
  * the optimal n is the eigenvector corresponding to the largest eigenvalue
- * of the matrix R*R^T, where R is the matrix whose i-th column is
- * r[idx[i]] - c.
+ * of the matrix R*R^T, where the i-th column of R is r[idx[i]] - c.
  *
- * @param[in]       m       Number of selected points (length of idx).
  * @param[in]       grid    Coordinates of all grid points.
  *                          grid[3*j], grid[3*j+1], grid[3*j+2] are the
  *                          x, y, z coordinates of the j-th point.
- *                          The length of grid is at least 3*m.
  * @param[in,out]   idx     Indices of the selected points within grid.
- *                          On exit, the indices are rearranged such that
- *                          points in each subset are put together.
+ *                          On return, idx will be rearranged such that
+ *                          points belonging to the same subset have their
+ *                          indices placed together.
+ * @param[in]       m       Number of selected points (length of idx).
  *
  * @return The number of points in the first subset (according to idx).
  *
  */
-int _maxmin_divide(int m, const double* grid, int* idx) {
+int _maxmin_divide(const double* grid, int* idx, int m) {
     assert(m > 1);
     if (m == 2) {
         return 1;
@@ -84,16 +83,16 @@ int _maxmin_divide(int m, const double* grid, int* idx) {
 
     int *head = idx;
     std::reverse_iterator<int*> tail(idx + m), rend(idx);
-    auto is_negative = [&dist](int j) { return dist[j] < 0; };
+    auto is_negative = [&dist, &idx](int& j) { return dist[&j - idx] < 0; };
     while ( ( head = std::find_if(head, idx + m, is_negative) ) <
             ( tail = std::find_if_not(tail, rend, is_negative) ).base() ) {
         std::swap(*head, *tail);
-        std::swap(dist[head - idx], dist[tail.base() - idx]);
+        std::swap(dist[head - idx], dist[tail.base() - idx - 1]);
         ++head;
         ++tail;
     }
 
-    return std::find(idx, idx + m, is_negative) - idx;
+    return head - idx;
 }
 
 } // end of anonymous namespace
@@ -109,7 +108,7 @@ std::vector<int> Grid::Batch::maxmin(
         return std::vector<int>{0};
     }
 
-    int m_left = _maxmin_divide(m, grid, idx);
+    int m_left = _maxmin_divide(grid, idx, m);
 
     std::vector<int> left = maxmin(grid, idx, m_left, m_thr);
     std::vector<int> right = maxmin(grid, idx + m_left, m - m_left, m_thr);

source/module_base/grid/batch.h

@@ -9,19 +9,39 @@ namespace Batch {
 /**
  * @brief Divide a set of points into batches by the "MaxMin" algorithm.
  *
- * This function recursively uses a cut plane to divide a set of grid
- * points into two subsets using the "MaxMin" algorithm, until the
- * number of points in each subset (batch) is no more than m_thr.
+ * This function recursively uses cut planes to divide grid points into
+ * two subsets using the "MaxMin" algorithm, until the number of points
+ * in each subset (batch) is no more than m_thr.
  *
  * @param[in]       grid        Coordinates of all grid points.
+ *                              grid[3*j], grid[3*j+1], grid[3*j+2] are
+ *                              the x, y, z coordinates of the j-th point.
  * @param[in,out]   idx         Indices of the initial set within grid.
  *                              On return, idx will be rearranged such
- *                              that points belonging to the same subset
- *                              are grouped together.
+ *                              that points belonging to the same batch
+ *                              have their indices placed together.
  * @param[in]       m           Number of points in the initial set.
- * @param[in]       m_thr       Size limit of subset.
+ *                              (length of idx)
+ * @param[in]       m_thr       Size limit of a batch.
  *
- * @return          Indices (within idx) of the first point in each batch.
+ * @return          Indices (for idx) of the first point in each batch.
+ *
+ * For example, given grid (illustrated by their indices) located as follows:
+ *
+ *      0  1  16          2  3  18
+ *      4  5  17            6  7
+ *
+ *
+ *      8  9             20 10 11
+ *     12 13 19           14 15
+ *
+ * a possible outcome with m_thr = 4 and idx(in) = {0, 1, 2, ..., 15} is:
+ *
+ * idx(out): 0, 1, 4, 5, 8, 9, 12, 13, 2, 3, 6, 7, 10, 11, 14, 15
+ * return  : {0, 4, 8, 12}
+ *
+ * which means the selected set (labeled 0-15) is divided into 4 batches:
+ * {0, 1, 4, 5}, {8, 9, 12, 13}, {2, 3, 6, 7}, {10, 11, 14, 15}.
  *
  */
 std::vector<int> maxmin(const double* grid, int* idx, int m, int m_thr);

source/module_base/grid/test/CMakeLists.txt

@@ -20,3 +20,10 @@ AddTest(
   ../radial.cpp
   ../delley.cpp
 )
+
+AddTest(
+  TARGET test_batch
+  SOURCES test_batch.cpp
+  ../batch.cpp
+  LIBS ${math_libs}
+)

source/module_base/grid/test/test_batch.cpp

@@ -10,7 +10,6 @@
 
 using namespace Grid::Batch;
 
-
 class BatchTest: public ::testing::Test
 {
 protected:
@@ -19,17 +18,28 @@ class BatchTest: public ::testing::Test
     std::vector<double> grid_;
     std::vector<int> idx_;
 
+    // parameters for cluster generation
     int n_each_ = 10;
-    double offset_ = 10.0;
     double width_ = 1.0;
+
+    // These offsets should be different from each other as maxmin might
+    // fail for highly symmetric, well-separated clusters.
+    // Consider the case where the 8 clusters as a whole have octahedral
+    // symmetry. In this case, R*R^T must be proprotional to the identity,
+    // and eigenvalues are three-fold degenerate, because xy, yz and zx
+    // plane are equivalent in terms of the maxmin optimization problem.
+    // This means eigenvectors are arbitrary in this case.
+    double offset_x_ = 7.0;
+    double offset_y_ = 8.0;
+    double offset_z_ = 9.0;
 };
 
-std::vector<double> gen_octant_cluster(int n_each, double offset, double width) {
+std::vector<double> gen_octant_cluster(int n_each, double offset_x, double offset_y, double offset_z, double width) {
 
     // Generates a set of points consisting of 8 well-separated, equal-sized
     // clusters located in individual octants.
 
-    std::vector<double> grid(n_each * 8);
+    std::vector<double> grid(n_each * 8 * 3);
     int I = 0;
 
     std::random_device rd;
@@ -40,15 +50,14 @@ std::vector<double> gen_octant_cluster(int n_each, double offset, double width)
         for (int sign_y : {-1, 1}) {
             for (int sign_z : {-1, 1}) {
                 for (int i = 0; i < n_each; ++i) {
-                    grid[3*I    ] = sign_x * offset + dis(gen);
-                    grid[3*I + 1] = sign_y * offset + dis(gen);
-                    grid[3*I + 2] = sign_z * offset + dis(gen);
+                    grid[3*I    ] = sign_x * offset_x + dis(gen);
+                    grid[3*I + 1] = sign_y * offset_y + dis(gen);
+                    grid[3*I + 2] = sign_z * offset_z + dis(gen);
                     ++I;
                 }
             }
         }
     }
-
     return grid;
 }
 
@@ -73,9 +82,9 @@ bool is_same_octant(int ngrid, const double* grid) {
 
 void BatchTest::SetUp()
 {
-    grid_ = gen_octant_cluster(n_each_, offset_, width_);
+    grid_ = gen_octant_cluster(n_each_, offset_x_, offset_y_, offset_z_, width_);
 
-    idx_.resize(grid_.size());
+    idx_.resize(grid_.size() / 3);
     std::iota(idx_.begin(), idx_.end(), 0);
 
     std::random_device rd;
@@ -91,21 +100,25 @@ TEST_F(BatchTest, MaxMinOctantCluster)
     // The resulting batches should be able to recover this structure.
 
     std::vector<int> delim = 
-        maxmin(n_each_, grid_.size(), grid_.data(), idx_.data());
+        maxmin(grid_.data(), idx_.data(), grid_.size() / 3, n_each_);
+
+    EXPECT_EQ(delim.size(), 8);
 
-    EXPECT_EQ(delim.size(), 7);
-    for (int i = 0; i < 7; ++i) {
-        // check number of points in each batch via index delimiters
-        EXPECT_EQ(delim[i], (i+1) * n_each_);
+    std::vector<double> grid_batch(3 * n_each_);
+    for (int i = 0; i < 8; ++i) {
 
-        // verify that points in each batch is in the same octant
-        std::vector<double> batch(3 * n_each_);
+        EXPECT_EQ(delim[i], i * n_each_);
+
+        // collect points within the present batch
         for (int j = 0; j < n_each_; ++j) {
-            for (int k = 0; k < 3; ++k) {
-                batch[3*j + k] = grid_[3*(i*n_each_ + j) + k];
-            }
+            int ig = idx_[delim[i] + j];
+            grid_batch[3*j    ] = grid_[3*ig    ];
+            grid_batch[3*j + 1] = grid_[3*ig + 1];
+            grid_batch[3*j + 2] = grid_[3*ig + 2];
         }
-        EXPECT_TRUE(is_same_octant(n_each_, batch.data()));
+
+        // verify that points in a batch reside in the same octant
+        EXPECT_TRUE(is_same_octant(n_each_, grid_batch.data()));
     }
 }