diff --git a/src/algorithm/nsga3.cpp b/src/algorithm/nsga3.cpp index 2fe54827..d1f35e03 100644 --- a/src/algorithm/nsga3.cpp +++ b/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 #include #include @@ -98,7 +97,7 @@ namespace gapp::algorithm }; RefLineGenerator ref_generator_; - detail::ConeTree ref_lines_; + std::vector ref_lines_; std::vector sol_info_; std::vector 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 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()); diff --git a/src/utility/algorithm.hpp b/src/utility/algorithm.hpp index dc4d6212..23ee91c7 100644 --- a/src/utility/algorithm.hpp +++ b/src/utility/algorithm.hpp @@ -144,6 +144,50 @@ namespace gapp::detail } } + template + requires std::invocable> + 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(value); + max_elem = first; + } + + return max_elem; + } + + template + requires std::invocable> + 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(value); + min_elem = first; + } + + return min_elem; + } + template constexpr void partial_shuffle(Iter first, Iter middle, Iter last, URBG&& gen) { diff --git a/src/utility/cone_tree.cpp b/src/utility/cone_tree.cpp index 756c6ba3..76d05f75 100644 --- a/src/utility/cone_tree.cpp +++ b/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; - - 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; - - 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(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)) { diff --git a/src/utility/cone_tree.hpp b/src/utility/cone_tree.hpp index 66a54c6e..a3cefd72 100644 --- a/src/utility/cone_tree.hpp +++ b/src/utility/cone_tree.hpp @@ -63,7 +63,7 @@ namespace gapp::detail Matrix points_; std::vector 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(); diff --git a/test/benchmark/cone_tree.cpp b/test/benchmark/cone_tree.cpp new file mode 100644 index 00000000..d6b3ebe6 --- /dev/null +++ b/test/benchmark/cone_tree.cpp @@ -0,0 +1,87 @@ +/* Copyright (c) 2023 Krisztián Rugási. Subject to the MIT License. */ + +#include +#include +#include +#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 randomPoints(size_t n, size_t dim) +{ + std::vector points(n); + for (Point& point : points) point = randomPoint(dim); + return points; +} + +static auto linearFind(const std::vector& 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)); }; +}