add Partition_N

Author: PierreMarchand20

CHANGELOG.md

@@ -48,18 +48,18 @@ All notable changes to this project will be documented in this file.
   - `task_based_internal_add_hmatrix_hmatrix_product` for task based alternative to `{sequential,openmp}_internal_add_hmatrix_hmatrix_product`.
   - `task_based_internal_triangular_hmatrix_hmatrix_solve` for task based alternative to `internal_triangular_hmatrix_hmatrix_solve`.
   - `task_based_lu_factorization` and `task_based_cholesky_factorization` for task based alternatives to `{sequential,openmp}_lu_factorization` and `{sequential,openmp}_cholesky_factorization`.
-  - `test_task_based_hmatrix_***.hpp` for testing various task based features.
+  - `test_task_based_hmatrix_*.hpp` for testing various task based features.
 - `internal_add_lrmat_hmatrix` is now overloaded to handle the case where the HMatrix is larger than the LowRankMatrix.
 - `get_leaves_from` is overloaded to return non const arguments.
 - `get_false_positive` in a tree builder.
 - `left_hmatrix_ancestor_of_right_hmatrix` and `left_hmatrix_descendant_of_right_hmatrix` for returning parent and children of a hmatrix.
+- `Partition_N` is an alternative to `Partition` for defining the partition of a cluster. The latter only splits along the principal axis of the cluster, while the former tries to be smarter.
 
 ### Changed
 
 - `VirtualInternalLowRankGenerator` and `VirtualLowRankGenerator`'s `copy_low_rank_approximation` function takes a `LowRankMatrix` as input to populate it and returns a boolean. The return value is true if the compression succeded, false otherwise.
 - `LowRankMatrix` constructors changed. It only takes sizes and an epsilon or a required rank. Then, it is expected to call a `VirtualInternalLowRankGenerator` to populate it.
 - `ClusterTreeBuilder` has now one strategy as `VirtualPartitioning`. Usual implementations are still available, for example using `Partitioning<double,ComputeLargestExtent,RegularSplitting>`.
-- When using `ClusterTreeBuilder` with `number_of_children=2^spatial_dimension`, it will do a binary/quad/octo-tree instead of `number_of_children` cut along the main direction.
 - `ClusterTreeBuilder` parameter `minclustersize` was removed, and a parameter `maximal_leaf_size` has been added.
 - `DistributedOperator` supports now both "global-to-local" and "local-to-local" operators, using respectively `VirtualGlobalToLocalOperator` and `VirtualLocalToLocalOperator` interfaces. The linear algebra associated has been updated to follow a more Blas-like interface.
 - `MatrixView` has been added to ease the use of matrix product. Most public functions for matrix products are also new templated to accept, `Matrix`, `MatrixView` or any other type following the same interface.

include/htool/clustering/implementations/partitioning.hpp

@@ -20,48 +20,129 @@ class Partitioning : public VirtualPartitioning<CoordinatePrecision> {
 
         // Compute directions
         Matrix<CoordinatePrecision> directions;
-        directions = ComputationDirectionPolicy::compute_direction(current_cluster, spatial_dimension, coordinates, radii, weights);
-
-        // If number of partitions corresponds to 2^d, we do binary/quadtree/octree/...
-        if (number_of_partitions == pow(2, spatial_dimension)) {
-            std::vector<int> dimensions(spatial_dimension);
-            std::iota(dimensions.begin(), dimensions.end(), int(0));
-            std::stack<std::tuple<int, int, int>> stack;
-            stack.push(std::make_tuple(current_cluster.get_offset(), current_cluster.get_size(), 0));
-            std::vector<std::pair<int, int>> result;
-            while (!stack.empty()) {
-                auto tmp_partition = stack.top();
-                stack.pop();
-                auto tmp_offset    = std::get<0>(tmp_partition);
-                auto tmp_size      = std::get<1>(tmp_partition);
-                auto tmp_dimension = std::get<2>(tmp_partition);
-
-                std::vector<CoordinatePrecision> direction = get_col(directions, tmp_dimension);
-
-                std::sort(permutation.begin() + tmp_offset, permutation.begin() + tmp_offset + tmp_size, [&](int a, int b) {
-                    CoordinatePrecision c = std::inner_product(coordinates + spatial_dimension * a, coordinates + spatial_dimension * (1 + a), direction.data(), CoordinatePrecision(0));
-                    CoordinatePrecision d = std::inner_product(coordinates + spatial_dimension * b, coordinates + spatial_dimension * (1 + b), direction.data(), CoordinatePrecision(0));
-                    return c < d;
-                });
+        std::vector<CoordinatePrecision> directions_weights;
+        std::tie(directions, directions_weights) = ComputationDirectionPolicy::compute_direction(current_cluster, spatial_dimension, coordinates, radii, weights);
 
-                auto tmp = SplittingPolicy::splitting(tmp_offset, tmp_size, spatial_dimension, coordinates, current_cluster.get_permutation(), direction, 2);
+        // Sort along main direction
+        std::sort(permutation.begin() + current_offset, permutation.begin() + current_offset + current_size, [&](int a, int b) {
+            CoordinatePrecision c = std::inner_product(coordinates + spatial_dimension * a, coordinates + spatial_dimension * (1 + a), directions.data(), CoordinatePrecision(0));
+            CoordinatePrecision d = std::inner_product(coordinates + spatial_dimension * b, coordinates + spatial_dimension * (1 + b), directions.data(), CoordinatePrecision(0));
+            return c < d;
+        });
 
-                if ((tmp_dimension < spatial_dimension - 1) and tmp.size() == 2) {
-                    stack.push(std::make_tuple(tmp[1].first, tmp[1].second, tmp_dimension + 1));
-                    stack.push(std::make_tuple(tmp[0].first, tmp[0].second, tmp_dimension + 1));
-                } else if ((tmp_dimension == spatial_dimension - 1) and tmp.size() == 2) {
-                    result.insert(result.end(), tmp.begin(), tmp.end());
-                } else {
-                    break;
-                }
+        // Split along main direction
+        return SplittingPolicy::splitting(current_cluster.get_offset(), current_cluster.get_size(), spatial_dimension, coordinates, current_cluster.get_permutation(), get_col(directions, 0), number_of_partitions);
+    }
+};
+
+template <typename CoordinatePrecision, typename ComputationDirectionPolicy, typename SplittingPolicy>
+class Partitioning_N : public VirtualPartitioning<CoordinatePrecision> {
+
+    // Generate all ordered remaining_d-decomposition of remaining_n
+    void backtrack(int remaining_n, int remaining_d, int start, std::vector<int> &current, std::vector<std::vector<int>> &results) {
+        if (remaining_d == 1) {
+            if (remaining_n <= start && remaining_n >= 1) {
+                current.push_back(remaining_n);
+                results.push_back(current);
+                current.pop_back();
+            }
+            return;
+        }
+
+        for (int f = start; f >= 1; f--) {
+            if (remaining_n % f == 0) {
+                current.push_back(f);
+                backtrack(remaining_n / f, remaining_d - 1, f, current, results);
+                current.pop_back();
+            }
+        }
+    }
+
+    CoordinatePrecision aspect_ratio(std::vector<CoordinatePrecision>) {
+    }
+
+    std::vector<int> distributed_splittings(int number_of_dimension, int number_of_partitions, std::vector<CoordinatePrecision> weights) {
+        // Generate all ordered integer decompositions of number_of_dimension elements
+        std::vector<std::vector<int>> integer_decompositions;
+        std::vector<int> current;
+        backtrack(number_of_partitions, number_of_dimension, number_of_partitions, current, integer_decompositions);
+
+        // Pick best decomposition for aspect ratio
+        CoordinatePrecision cost = std::numeric_limits<CoordinatePrecision>::max();
+        int index                = 0;
+        int result_index;
+        std::vector<CoordinatePrecision> aspect_ratios(number_of_dimension);
+        for (auto &integer_decomposition : integer_decompositions) {
+            std::transform(integer_decomposition.cbegin(), integer_decomposition.cend(), weights.cbegin(), aspect_ratios.begin(), [](const CoordinatePrecision &p, const CoordinatePrecision &a) { return a / p; });
+            CoordinatePrecision current_cost = *std::max_element(aspect_ratios.begin(), aspect_ratios.end()) / *std::min_element(aspect_ratios.begin(), aspect_ratios.end());
+
+            if (current_cost < cost) {
+                cost         = current_cost;
+                result_index = index;
+            }
+            index++;
+        }
+        return integer_decompositions[result_index];
+    }
+
+  public:
+    std::vector<std::pair<int, int>> compute_partitioning(Cluster<CoordinatePrecision> &current_cluster, int spatial_dimension, const CoordinatePrecision *coordinates, const CoordinatePrecision *const radii, const CoordinatePrecision *const weights, int number_of_partitions) override {
+
+        std::vector<int> &permutation = current_cluster.get_permutation();
+        auto current_offset           = current_cluster.get_offset();
+        auto current_size             = current_cluster.get_size();
+
+        // Compute directions
+        Matrix<CoordinatePrecision> directions;
+        std::vector<CoordinatePrecision> directions_weights;
+        std::tie(directions, directions_weights) = ComputationDirectionPolicy::compute_direction(current_cluster, spatial_dimension, coordinates, radii, weights);
+
+        int number_of_relevant_directions = 0;
+        for (auto direction_weight : directions_weights) {
+            if (direction_weight > std::numeric_limits<CoordinatePrecision>::epsilon() * 10) {
+                number_of_relevant_directions += 1;
             }
+        }
+
+        std::vector<int> number_of_splittings_by_direction = distributed_splittings(number_of_relevant_directions, number_of_partitions, directions_weights);
+        number_of_relevant_directions                      = number_of_splittings_by_direction.size();
+
+        std::stack<std::tuple<int, int, int>> stack;
+        stack.push(std::make_tuple(current_cluster.get_offset(), current_cluster.get_size(), 0));
+        std::vector<std::pair<int, int>> result;
+        while (!stack.empty()) {
+            auto tmp_partition = stack.top();
+            stack.pop();
+            auto tmp_offset    = std::get<0>(tmp_partition);
+            auto tmp_size      = std::get<1>(tmp_partition);
+            auto tmp_dimension = std::get<2>(tmp_partition);
+
+            std::vector<CoordinatePrecision> direction = get_col(directions, tmp_dimension);
 
-            if (result.size() == number_of_partitions) {
-                std::sort(result.begin(), result.end(), [](auto a, auto b) { return a.first < b.first; });
-                return result;
+            std::sort(permutation.begin() + tmp_offset, permutation.begin() + tmp_offset + tmp_size, [&](int a, int b) {
+                CoordinatePrecision c = std::inner_product(coordinates + spatial_dimension * a, coordinates + spatial_dimension * (1 + a), direction.data(), CoordinatePrecision(0));
+                CoordinatePrecision d = std::inner_product(coordinates + spatial_dimension * b, coordinates + spatial_dimension * (1 + b), direction.data(), CoordinatePrecision(0));
+                return c < d;
+            });
+
+            auto tmp = SplittingPolicy::splitting(tmp_offset, tmp_size, spatial_dimension, coordinates, current_cluster.get_permutation(), direction, number_of_splittings_by_direction[tmp_dimension]);
+
+            if ((tmp_dimension < number_of_relevant_directions - 1) and tmp.size() == number_of_splittings_by_direction[tmp_dimension]) {
+                for (int p = number_of_splittings_by_direction[tmp_dimension] - 1; p >= 0; p--) {
+                    stack.push(std::make_tuple(tmp[p].first, tmp[p].second, tmp_dimension + 1));
+                }
+            } else if ((tmp_dimension == number_of_relevant_directions - 1) and tmp.size() == number_of_splittings_by_direction[tmp_dimension]) {
+                result.insert(result.end(), tmp.begin(), tmp.end());
+            } else {
+                break;
             }
         }
 
+        if (result.size() == number_of_partitions) {
+            std::sort(result.begin(), result.end(), [](auto a, auto b) { return a.first < b.first; });
+            return result;
+        }
+
         // Sort along main direction
         std::sort(permutation.begin() + current_offset, permutation.begin() + current_offset + current_size, [&](int a, int b) {
             CoordinatePrecision c = std::inner_product(coordinates + spatial_dimension * a, coordinates + spatial_dimension * (1 + a), directions.data(), CoordinatePrecision(0));
@@ -77,14 +158,13 @@ class Partitioning : public VirtualPartitioning<CoordinatePrecision> {
 template <typename T>
 class ComputeLargestExtent {
   public:
-    static Matrix<T> compute_direction(const Cluster<T> &cluster, int spatial_dimension, const T *const coordinates, const T *const, const T *const weights) {
+    static std::pair<Matrix<T>, std::vector<T>> compute_direction(const Cluster<T> &cluster, int spatial_dimension, const T *const coordinates, const T *const, const T *const weights) {
         if (spatial_dimension != 2 && spatial_dimension != 3) {
             htool::Logger::get_instance().log(LogLevel::ERROR, "clustering not define for spatial dimension !=2 and !=3"); // LCOV_EXCL_LINE
         }
 
         const std::vector<int> &permutation = cluster.get_permutation();
         Matrix<T> cov(spatial_dimension, spatial_dimension);
-        Matrix<T> directions(spatial_dimension, spatial_dimension);
 
         for (int j = 0; j < cluster.get_size(); j++) {
             std::vector<T> u(spatial_dimension, 0);
@@ -98,19 +178,23 @@ class ComputeLargestExtent {
                 }
             }
         }
-        if (spatial_dimension == 2) {
-            directions = solve_EVP_2(cov);
-        } else if (spatial_dimension == 3) {
-            directions = solve_EVP_3(cov);
+
+        auto evp_eigs = spatial_dimension == 2 ? solve_EVP_2(cov) : solve_EVP_3(cov);
+        for (auto &eig : evp_eigs.second) {
+            if (eig > 0) {
+                eig = std::sqrt(eig);
+            } else {
+                eig = 0;
+            }
         }
-        return directions;
+        return evp_eigs;
     }
 };
 
 template <typename T>
 class ComputeBoundingBox {
   public:
-    static Matrix<T> compute_direction(const Cluster<T> &cluster, int spatial_dimension, const T *const coordinates, const T *const, const T *const) {
+    static std::pair<Matrix<T>, std::vector<T>> compute_direction(const Cluster<T> &cluster, int spatial_dimension, const T *const coordinates, const T *const, const T *const) {
 
         const std::vector<int> &permutation = cluster.get_permutation();
 
@@ -131,14 +215,17 @@ class ComputeBoundingBox {
         }
 
         // Direction of largest extent
+        std::vector<T> lengths(spatial_dimension);
         std::vector<int> indexes(spatial_dimension);
         std::iota(indexes.begin(), indexes.end(), int(0));
         std::sort(indexes.begin(), indexes.end(), [&max_point, &min_point](int a, int b) { return (max_point[a] - min_point[a]) < (max_point[b] - min_point[b]); });
         Matrix<T> directions(spatial_dimension, spatial_dimension);
         for (int dim = 0; dim < spatial_dimension; dim++) {
             directions(indexes[spatial_dimension - 1 - dim], dim) = 1;
+
+            lengths[dim] = max_point[indexes[spatial_dimension - 1 - dim]] - min_point[indexes[spatial_dimension - 1 - dim]];
         }
-        return directions;
+        return std::make_pair(directions, lengths);
     }
 };
 

include/htool/misc/evp.hpp

@@ -10,7 +10,7 @@
 namespace htool {
 
 template <typename T>
-Matrix<T> solve_EVP_2(const Matrix<T> &cov) {
+std::pair<Matrix<T>, std::vector<T>> solve_EVP_2(const Matrix<T> &cov) {
     std::vector<T> dir(2, 0);
     std::vector<T> eigs(2);
     Matrix<T> I(2, 2);
@@ -45,11 +45,11 @@ Matrix<T> solve_EVP_2(const Matrix<T> &cov) {
         result(0, 0) = 1;
         result(1, 1) = 1;
     }
-    return result;
+    return std::make_pair(result, eigs);
 }
 
 template <typename T>
-Matrix<T> solve_EVP_3(const Matrix<T> &cov) {
+std::pair<Matrix<T>, std::vector<T>> solve_EVP_3(const Matrix<T> &cov) {
     Matrix<T> result(3, 3);
     T p1 = std::pow(cov(0, 1), 2) + std::pow(cov(0, 2), 2) + std::pow(cov(1, 2), 2);
     std::vector<T> eigs(3);
@@ -150,7 +150,7 @@ Matrix<T> solve_EVP_3(const Matrix<T> &cov) {
             result(2, 2) = 1;
         }
     }
-    return result;
+    return std::make_pair(result, eigs);
 }
 } // namespace htool
 #endif

include/htool/testing/geometry.hpp

@@ -8,25 +8,40 @@
 namespace htool {
 
 template <typename T>
-void create_disk(int space_dim, T z, int nr, T *const xt) {
-
+void create_rotated_ellipse(int space_dim, T a, T b, T alpha, T z, int nr, T *const xt) {
     std::mt19937 mersenne_engine(0);
     std::uniform_real_distribution<T> dist(0, 1);
     auto gen = [&dist, &mersenne_engine]() {
         return dist(mersenne_engine);
     };
-    T z1 = z;
+
+    T cos_alpha = std::cos(alpha);
+    T sin_alpha = std::sin(alpha);
+
     for (int j = 0; j < nr; j++) {
-        T rho                 = gen(); // (T) otherwise integer division!
-        T theta               = gen();
-        xt[space_dim * j + 0] = std::sqrt(rho) * std::cos(2 * static_cast<T>(M_PI) * theta);
-        xt[space_dim * j + 1] = std::sqrt(rho) * std::sin(2 * static_cast<T>(M_PI) * theta);
+        T rho   = gen();
+        T theta = gen();
+        T r     = std::sqrt(rho);
+        T phi   = 2 * static_cast<T>(M_PI) * theta;
+
+        // Axis-aligned ellipse coordinates
+        T x_prime = a * r * std::cos(phi);
+        T y_prime = b * r * std::sin(phi);
+
+        // Apply rotation
+        xt[space_dim * j + 0] = cos_alpha * x_prime - sin_alpha * y_prime;
+        xt[space_dim * j + 1] = sin_alpha * x_prime + cos_alpha * y_prime;
+
         if (space_dim == 3)
-            xt[space_dim * j + 2] = z1;
-        // sqrt(rho) otherwise the points would be concentrated in the center of the disk
+            xt[space_dim * j + 2] = z;
     }
 }
 
+template <typename T>
+void create_disk(int space_dim, T z, int nr, T *const xt) {
+    create_rotated_ellipse(space_dim, T(1.), T(1.), T(0.), z, nr, xt);
+}
+
 template <typename T>
 void create_sphere(int nr, T *const xt, std::array<T, 3> offset = {0, 0, 0}) {
 

include/htool/testing/partition.hpp

@@ -14,6 +14,7 @@ template <typename CoordinatePrecision>
 std::vector<int> test_global_partition(int spatial_dimension, int number_of_points, const std::vector<CoordinatePrecision> &coordinates, int partition_size) {
     // Compute largest extent
     Matrix<CoordinatePrecision> direction(spatial_dimension, spatial_dimension);
+    std::vector<CoordinatePrecision> weigths;
     Matrix<CoordinatePrecision> cov(spatial_dimension, spatial_dimension);
     std::vector<int> permutation(number_of_points, 0);
     std::iota(permutation.begin(), permutation.end(), int(0));
@@ -31,9 +32,9 @@ std::vector<int> test_global_partition(int spatial_dimension, int number_of_poin
         }
     }
     if (spatial_dimension == 2) {
-        direction = solve_EVP_2(cov);
+        std::tie(direction, weigths) = solve_EVP_2(cov);
     } else if (spatial_dimension == 3) {
-        direction = solve_EVP_3(cov);
+        std::tie(direction, weigths) = solve_EVP_3(cov);
     }
 
     // sort
@@ -62,6 +63,7 @@ template <typename CoordinatePrecision>
 std::vector<int> test_local_partition(int spatial_dimension, int number_of_points, std::vector<CoordinatePrecision> &coordinates, int partition_size) {
     // Compute largest extent
     Matrix<CoordinatePrecision> direction(spatial_dimension, spatial_dimension);
+    std::vector<CoordinatePrecision> weigths;
     Matrix<CoordinatePrecision> cov(spatial_dimension, spatial_dimension);
     std::vector<int> permutation(number_of_points, 0);
     std::iota(permutation.begin(), permutation.end(), int(0));
@@ -79,9 +81,9 @@ std::vector<int> test_local_partition(int spatial_dimension, int number_of_point
         }
     }
     if (spatial_dimension == 2) {
-        direction = solve_EVP_2(cov);
+        std::tie(direction, weigths) = solve_EVP_2(cov);
     } else if (spatial_dimension == 3) {
-        direction = solve_EVP_3(cov);
+        std::tie(direction, weigths) = solve_EVP_3(cov);
     }
 
     // sort