add ConeTree benchmarks & use linear search in the NSGA3

src/algorithm/nsga3.cpp

@@ -11,7 +11,6 @@
 #include "../utility/math.hpp"
 #include "../utility/rng.hpp"
 #include "../utility/utility.hpp"
-#include "../utility/cone_tree.hpp"
 #include <algorithm>
 #include <execution>
 #include <functional>
@@ -98,7 +97,7 @@ namespace gapp::algorithm
         };
 
         RefLineGenerator ref_generator_;
-        detail::ConeTree ref_lines_;
+        std::vector<Point> ref_lines_;
 
         std::vector<CandidateInfo> sol_info_;
         std::vector<size_t> niche_counts_;
@@ -200,7 +199,7 @@ namespace gapp::algorithm
         for (size_t dim = 0; dim < ideal_point_.size(); dim++)
         {
             auto weights = weightVector(ideal_point_.size(), dim);
-            auto ASFi = [&](const auto& fvec) noexcept { return ASF(ideal_point_, weights, fvec); };
+            auto ASFi = [&](const auto& fvec) { return ASF(ideal_point_, weights, fvec); };
 
             std::vector<double> chebysev_distances(popsize + extreme_points_.size());
 
@@ -244,10 +243,12 @@ namespace gapp::algorithm
         {
             const FitnessVector fnorm = normalizeFitnessVec(first[sol.idx], ideal_point_, nadir_point_);
 
-            const auto best = ref_lines_.findBestMatch(fnorm);
+            auto idistance = [&](const auto& line) { return std::inner_product(fnorm.begin(), fnorm.end(), line.begin(), 0.0); };
 
-            sol_info_[sol.idx].ref_idx = size_t(best.elem - ref_lines_.begin());
-            sol_info_[sol.idx].ref_dist = math::perpendicularDistanceSq(*best.elem, fnorm);
+            auto closest = detail::max_element(ref_lines_.begin(), ref_lines_.end(), idistance);
+
+            sol_info_[sol.idx].ref_idx  = std::distance(ref_lines_.begin(), closest);
+            sol_info_[sol.idx].ref_dist = math::perpendicularDistanceSq(*closest, fnorm);
         });
     }
 
@@ -349,8 +350,7 @@ namespace gapp::algorithm
         pimpl_->ideal_point_ = detail::maxFitness(fitness_matrix.begin(), fitness_matrix.end());
         pimpl_->extreme_points_ = {};
 
-        auto ref_lines = pimpl_->generateReferencePoints(ga.num_objectives(), ga.population_size());
-        pimpl_->ref_lines_ = detail::ConeTree{ ref_lines };
+        pimpl_->ref_lines_ = pimpl_->generateReferencePoints(ga.num_objectives(), ga.population_size());
         pimpl_->niche_counts_.resize(pimpl_->ref_lines_.size());
 
         auto pfronts = nonDominatedSort(fitness_matrix.begin(), fitness_matrix.end());

src/utility/algorithm.hpp

@@ -144,6 +144,50 @@ namespace gapp::detail
         }
     }
 
+    template<std::forward_iterator Iter, typename F>
+    requires std::invocable<F&, std::iter_reference_t<Iter>>
+    constexpr Iter max_element(Iter first, Iter last, F&& transform)
+    {
+        GAPP_ASSERT(std::distance(first, last) > 0);
+
+        if (first == last) return first;
+
+        Iter max_elem = first;
+        auto&& max_value = std::invoke(transform, *first);
+
+        while (++first != last)
+        {
+            auto&& value = std::invoke(transform, *first);
+            if (value <= max_value) continue;
+            max_value = std::forward<decltype(value)>(value);
+            max_elem = first;
+        }
+
+        return max_elem;
+    }
+
+    template<std::forward_iterator Iter, typename F>
+    requires std::invocable<F&, std::iter_reference_t<Iter>>
+    constexpr Iter min_element(Iter first, Iter last, F&& transform)
+    {
+        GAPP_ASSERT(std::distance(first, last) > 0);
+
+        if (first == last) return first;
+
+        Iter min_elem = first;
+        auto&& min_value = std::invoke(transform, *first);
+
+        while (++first != last)
+        {
+            auto&& value = std::invoke(transform, *first);
+            if (value >= min_value) continue;
+            min_value = std::forward<decltype(value)>(value);
+            min_elem = first;
+        }
+
+        return min_elem;
+    }
+
     template<std::random_access_iterator Iter, typename URBG>
     constexpr void partial_shuffle(Iter first, Iter middle, Iter last, URBG&& gen)
     {

src/utility/cone_tree.cpp

@@ -1,6 +1,7 @@
 /* Copyright (c) 2022 Krisztián Rugási. Subject to the MIT License. */
 
 #include "cone_tree.hpp"
+#include "algorithm.hpp"
 #include "functional.hpp"
 #include "math.hpp"
 #include "matrix.hpp"
@@ -40,23 +41,7 @@ namespace gapp::detail
     /* Find the point in the range [first, last) furthest from a point (using Euclidean distances). */
     static inline const_iterator findFurthestElement(const_iterator first, const_iterator last, const_iterator from)
     {
-        GAPP_ASSERT(std::distance(first, last) > 0);
-
-        const_iterator furthest;
-        double max_distance = -math::inf<double>;
-
-        for (; first != last; ++first)
-        {
-            const double distance = math::euclideanDistanceSq(*first, *from);
-
-            if (distance > max_distance)
-            {
-                furthest = first;
-                max_distance = distance;
-            }
-        }
-
-        return furthest;
+        return detail::max_element(first, last, std::bind_front(math::euclideanDistanceSq, *from));
     }
 
     /* Find the 2 partition points that will be used to split the range [first, last) into 2 parts. */
@@ -94,23 +79,16 @@ namespace gapp::detail
     /* Find the Euclidean distance between the center point and the point in the range [first, last) furthest from it. */
     static inline double findRadius(const_iterator first, const_iterator last, const Point& center)
     {
-        GAPP_ASSERT(std::distance(first, last) > 0);
-
-        double max_distance = -math::inf<double>;
-
-        for (; first != last; ++first)
-        {
-            const double distance = math::euclideanDistanceSq(*first, center);
-            max_distance = std::max(max_distance, distance);
-        }
+        auto distance = std::bind_front(math::euclideanDistanceSq, center);
+        auto furthest = detail::max_element(first, last, distance);
 
-        return std::sqrt(max_distance);
+        return std::sqrt(distance(*furthest));
     }
 
     /* Returns true if the node is a leaf node. */
     static inline bool isLeafNode(const Node& node) noexcept
     {
-        return !static_cast<bool>(node.left & node.right);
+        return (node.left == 0) && (node.right == 0);
     }
 
     /* Return the max possible inner product between the point and a point inside the node. */
@@ -223,6 +201,8 @@ namespace gapp::detail
             }
             else
             {
+                GAPP_ASSERT(node->left && node->right);
+
                 if (innerProductUpperBound(left_child(*node), query_point, query_norm) <
                     innerProductUpperBound(right_child(*node), query_point, query_norm))
                 {

src/utility/cone_tree.hpp

@@ -63,7 +63,7 @@ namespace gapp::detail
         Matrix<double> points_;
         std::vector<Node> nodes_;
 
-        inline static constexpr size_t MAX_LEAF_ELEMENTS = 22;  /* The maximum number of points in a leaf node. */
+        inline static constexpr size_t MAX_LEAF_ELEMENTS = 8;  /* The maximum number of points in a leaf node. */
 
         void buildTree();
 

test/benchmark/cone_tree.cpp

@@ -0,0 +1,87 @@
+/* Copyright (c) 2023 Krisztián Rugási. Subject to the MIT License. */
+
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/generators/catch_generators.hpp>
+#include <catch2/benchmark/catch_benchmark.hpp>
+#include "utility/cone_tree.hpp"
+#include "utility/rng.hpp"
+#include "utility/algorithm.hpp"
+
+using namespace gapp;
+using namespace gapp::detail;
+using namespace gapp::math;
+using namespace Catch;
+
+
+static Point randomPoint(size_t dim)
+{
+    Point point(dim);
+    for (double& elem : point) elem = rng::randomReal();
+    return point;
+}
+
+static std::vector<Point> randomPoints(size_t n, size_t dim)
+{
+    std::vector<Point> points(n);
+    for (Point& point : points) point = randomPoint(dim);
+    return points;
+}
+
+static auto linearFind(const std::vector<Point>& lines, const Point& point)
+{
+    auto idistance = [&](const auto& line) { return std::inner_product(point.begin(), point.end(), line.begin(), 0.0); };
+
+    return detail::max_element(lines.begin(), lines.end(), idistance);
+}
+
+
+TEST_CASE("cone_tree_ctor", "[cone_tree]")
+{
+    constexpr size_t ndim = 3;
+
+    BENCHMARK_ADVANCED("small")(Benchmark::Chronometer meter)
+    {
+        auto points = randomPoints(100, ndim);
+        meter.measure([&] { return ConeTree{ points }; });
+    };
+
+    BENCHMARK_ADVANCED("medium")(Benchmark::Chronometer meter)
+    {
+        auto points = randomPoints(1000, ndim);
+        meter.measure([&] { return ConeTree{ points }; });
+    };
+
+    BENCHMARK_ADVANCED("large")(Benchmark::Chronometer meter)
+    {
+        auto points = randomPoints(10000, ndim);
+        meter.measure([&] { return ConeTree{ points }; });
+    };
+}
+
+TEST_CASE("cone_tree_lookup_size", "[cone_tree]")
+{
+    constexpr size_t ndim = 3;
+    const size_t size = GENERATE(100, 1000, 10000);
+
+    WARN("Number of points: " << size);
+
+    auto points = randomPoints(size, ndim);
+    auto tree = ConeTree{ points };
+
+    BENCHMARK("cone_tree") { return tree.findBestMatch(randomPoint(ndim)); };
+    BENCHMARK("linsearch") { return linearFind(points, randomPoint(ndim)); };
+}
+
+TEST_CASE("cone_tree_lookup_dim", "[cone_tree]")
+{
+    const size_t ndim = GENERATE(3, 15, 100);
+    constexpr size_t size = 10000;
+
+    WARN("Number of dimensions: " << ndim);
+
+    auto points = randomPoints(size, ndim);
+    auto tree = ConeTree{ points };
+
+    BENCHMARK("cone_tree") { return tree.findBestMatch(randomPoint(ndim)); };
+    BENCHMARK("linsearch") { return linearFind(points, randomPoint(ndim)); };
+}