diff --git a/README.md b/README.md
index 5de06fff..a21de16d 100644
--- a/README.md
+++ b/README.md
@@ -74,7 +74,7 @@ will install all available configurations:
 
 ```shell
 git clone https://github.com/KRM7/gapp.git
-sudo gapp/build/install.sh
+sudo bash ./gapp/build/install.sh
 ```
 
 Once the library is installed, you can import it using `find_package` and then link
diff --git a/docs/metrics.md b/docs/metrics.md
index 86597b86..8ba1b315 100644
--- a/docs/metrics.md
+++ b/docs/metrics.md
@@ -64,20 +64,14 @@ are intended to be used for multi-objective optimization problems.
 
 If you want to track something that doesn't have a metric already implemented 
 for it, it's possible to implement your own metrics. Metrics must be derived
-from the `Monitor` class, and implement 3 methods: `initialize`, `update`, and
-`value_at`:
+from the `Monitor` class, and implement the `initialize`, and `update` methods:
 
 ```cpp
 // The second type parameter of Monitor is the type used to store the
 // gathered metrics. This will be used as the type of the data_ field.
 class MyMetric : public metrics::Monitor<MyMetric, std::vector<double>>
 {
-public:
-    // Returns the value of the metric in the given generation.
-    // Note that this method is not virtual.
-    double value_at(size_t generation) const noexcept { return data_[generation]; }
-private:
-    // Initialize the metric. Called at the start of a run.
+    // (optional) Initialize the metric. Called at the start of a run.
     void initialize(const GaInfo& ga) override { data_.clear(); }
 
     // Update the metric with a new value from the current generation.
diff --git a/examples/9_metrics.cpp b/examples/9_metrics.cpp
index cb72b928..b38d6775 100644
--- a/examples/9_metrics.cpp
+++ b/examples/9_metrics.cpp
@@ -10,8 +10,6 @@ using namespace gapp;
 
 struct MyMetric : public metrics::Monitor<MyMetric, std::vector<double>>
 {
-    double value_at(size_t generation) const noexcept { return data_[generation]; }
-    void initialize(const GaInfo&) override { data_.clear(); }
     void update(const GaInfo& ga) override { data_.push_back(ga.fitness_matrix()[0][0]); }
 };
 
@@ -29,6 +27,6 @@ int main()
         std::cout << std::format("Generation {}\t| {:.6f}\n", gen + 1, metric[gen]);
     }
 
-    const auto* hypervol = GA.get_metric_if<metrics::AutoHypervolume>(); // untracked metric
+    [[maybe_unused]] const auto* hypervol = GA.get_metric_if<metrics::AutoHypervolume>(); // untracked metric
     assert(hypervol == nullptr);
 }
diff --git a/src/algorithm/nd_sort.cpp b/src/algorithm/nd_sort.cpp
index e764fea7..3f23dea0 100644
--- a/src/algorithm/nd_sort.cpp
+++ b/src/algorithm/nd_sort.cpp
@@ -35,7 +35,7 @@ namespace gapp::algorithm::dtl
     {
         GAPP_ASSERT(std::distance(current, last) >= 0);
 
-        return std::find_if(current, last, compose(&FrontInfo::rank, detail::not_equal_to(current->rank)));
+        return std::ranges::find_if(current, last, detail::not_equal_to(current->rank), &FrontInfo::rank);
     }
 
     std::vector<ParetoFrontsRange> paretoFrontBounds(ParetoFronts& pareto_fronts)
@@ -65,8 +65,7 @@ namespace gapp::algorithm::dtl
         auto last_in = first + popsize;
 
         if (last_in == last) return { last, last };
-
-        auto partial_front_first = std::find_if(first, last_in, compose(&FrontInfo::rank, detail::equal_to(last_in->rank)));
+        auto partial_front_first = std::ranges::find_if(first, last_in, detail::equal_to(last_in->rank), &FrontInfo::rank);
         auto partial_front_last = nextFrontBegin(std::prev(last_in), last);
 
         return { partial_front_first, partial_front_last };
@@ -205,7 +204,10 @@ namespace gapp::algorithm::dtl
 
                 std::for_each(last_idx, indices.rend(), [&](size_t col) noexcept
                 {
-                    if (dmat(row, col).load(std::memory_order_relaxed) == MAXIMAL) dmat(row, col).store(NONMAXIMAL, std::memory_order_relaxed);
+                    if (dmat(row, col).load(std::memory_order_relaxed) == MAXIMAL)
+                    {
+                        dmat(row, col).store(NONMAXIMAL, std::memory_order_relaxed);
+                    }
                 });
             });
         });
diff --git a/src/algorithm/nsga3.cpp b/src/algorithm/nsga3.cpp
index ae01b35c..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 <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());
@@ -421,8 +421,7 @@ namespace gapp::algorithm
 
     std::vector<size_t> NSGA3::optimalSolutionsImpl(const GaInfo&) const
     {
-        return detail::find_indices(pimpl_->sol_info_,
-                                    detail::compose(&Impl::CandidateInfo::rank, detail::equal_to(0_sz)));
+        return detail::find_indices(pimpl_->sol_info_, [](const Impl::CandidateInfo& sol) { return sol.rank == 0; });
     }
 
 } // namespace gapp::algorithm
\ No newline at end of file
diff --git a/src/core/ga_base.impl.hpp b/src/core/ga_base.impl.hpp
index 31266af0..325c9b50 100644
--- a/src/core/ga_base.impl.hpp
+++ b/src/core/ga_base.impl.hpp
@@ -358,6 +358,8 @@ namespace gapp
 
         metrics_.initialize(*this);
         metrics_.update(*this);
+
+        if (end_of_generation_callback_) end_of_generation_callback_(*this);
     }
 
     template<typename T>
@@ -481,7 +483,16 @@ namespace gapp
         auto optimal_pop = algorithm_->optimalSolutions(*this, pop);
 
         optimal_sols = detail::mergeParetoSets(std::move(optimal_sols), std::move(optimal_pop));
-        detail::erase_duplicates(optimal_sols);
+
+        /* Duplicate elements are removed from optimal_sols using exact comparison
+        *  of the chromosomes in order to avoid issues with using a non-transitive
+        *  comparison function for std::sort and std::unique. */
+        auto chrom_eq = [](const auto& lhs, const auto& rhs) { return lhs.chromosome == rhs.chromosome; };
+        auto chrom_less = [](const auto& lhs, const auto& rhs) { return lhs.chromosome < rhs.chromosome; };
+
+        std::sort(optimal_sols.begin(), optimal_sols.end(), chrom_less);
+        auto last = std::unique(optimal_sols.begin(), optimal_sols.end(), chrom_eq);
+        optimal_sols.erase(last, optimal_sols.end());
     }
 
     template<typename T>
diff --git a/src/metrics/distribution_metrics.cpp b/src/metrics/distribution_metrics.cpp
index a8f8495a..07c72ab4 100644
--- a/src/metrics/distribution_metrics.cpp
+++ b/src/metrics/distribution_metrics.cpp
@@ -13,13 +13,6 @@ namespace gapp::metrics
 {
     using math::Point;
 
-    std::span<const double> NadirPoint::value_at(size_t generation) const noexcept
-    {
-        GAPP_ASSERT(generation < data_.size());
-
-        return data_[generation];
-    }
-
     void NadirPoint::initialize(const GaInfo& ga)
     {
         data_.clear();
@@ -39,13 +32,6 @@ namespace gapp::metrics
         ref_point_(std::move(ref_point))
     {}
 
-    double Hypervolume::value_at(size_t generation) const noexcept
-    {
-        GAPP_ASSERT(generation < data_.size());
-
-        return data_[generation];
-    }
-
     void Hypervolume::initialize(const GaInfo& ga)
     {
         GAPP_ASSERT(ref_point_.size() == ga.num_objectives());
@@ -63,13 +49,6 @@ namespace gapp::metrics
     }
 
 
-    double AutoHypervolume::value_at(size_t generation) const noexcept
-    {
-        GAPP_ASSERT(generation < data_.size());
-
-        return data_[generation];
-    }
-
     void AutoHypervolume::initialize(const GaInfo& ga)
     {
         data_.clear();
diff --git a/src/metrics/distribution_metrics.hpp b/src/metrics/distribution_metrics.hpp
index 312efb9b..7a51d36e 100644
--- a/src/metrics/distribution_metrics.hpp
+++ b/src/metrics/distribution_metrics.hpp
@@ -16,9 +16,6 @@ namespace gapp::metrics
     */
     class NadirPoint final : public Monitor<NadirPoint, FitnessMatrix>
     {
-    public:
-        std::span<const double> value_at(size_t generation) const noexcept;
-    private:
         void initialize(const GaInfo& ga) override;
         void update(const GaInfo& ga) override;
     };
@@ -50,8 +47,6 @@ namespace gapp::metrics
         /** @returns The reference point used for computing the hypervolumes. */
         const FitnessVector& ref_point() const noexcept { return ref_point_; }
 
-        double value_at(size_t generation) const noexcept;
-
     private:
         void initialize(const GaInfo& ga) override;
         void update(const GaInfo& ga) override;
@@ -79,7 +74,6 @@ namespace gapp::metrics
         /** @returns The reference point used for computing the hypervolumes. */
         const FitnessVector& ref_point() const noexcept { return worst_point_; }
 
-        double value_at(size_t generation) const noexcept;
     private:
         void initialize(const GaInfo& ga) override;
         void update(const GaInfo& ga) override;
diff --git a/src/metrics/fitness_metrics.hpp b/src/metrics/fitness_metrics.hpp
index a5a405e0..32de34ce 100644
--- a/src/metrics/fitness_metrics.hpp
+++ b/src/metrics/fitness_metrics.hpp
@@ -16,13 +16,6 @@ namespace gapp::metrics
     template<typename Derived>
     class FitnessMonitor : public Monitor<Derived, FitnessMatrix>
     {
-    public:
-        std::span<const double> value_at(size_t generation) const noexcept
-        {
-            GAPP_ASSERT(generation < this->data_.size());
-
-            return this->data_[generation];
-        }
     private:
         void initialize(const GaInfo& ga) override
         {
diff --git a/src/metrics/metric_set.hpp b/src/metrics/metric_set.hpp
index 075d75c2..9732737a 100644
--- a/src/metrics/metric_set.hpp
+++ b/src/metrics/metric_set.hpp
@@ -51,7 +51,7 @@ namespace gapp::detail
     MetricSet::MetricSet(Metrics... metrics)
     {
         metrics_.reserve(sizeof...(metrics));
-        (metrics_.push_back(std::make_unique<Metrics>(std::move(metrics))), ...);
+        ( metrics_.push_back(std::make_unique<Metrics>(std::move(metrics))), ... );
     }
 
     template<typename Metric>
@@ -60,8 +60,7 @@ namespace gapp::detail
     {
         auto found = std::find_if(metrics_.begin(), metrics_.end(), [](const auto& metric) { return metric->type_id() == detail::type_id<Metric>; });
 
-        if (found == metrics_.end()) return nullptr;
-        return static_cast<const Metric*>(found->get());
+        return found != metrics_.end() ? static_cast<Metric*>(found->get()) : nullptr;
     }
 
 } // namespace gapp::detail
diff --git a/src/metrics/misc_metrics.cpp b/src/metrics/misc_metrics.cpp
index 25f662fa..f3c6145f 100644
--- a/src/metrics/misc_metrics.cpp
+++ b/src/metrics/misc_metrics.cpp
@@ -22,11 +22,4 @@ namespace gapp::metrics
         data_.push_back(sum_ - old_sum);
     }
 
-    size_t FitnessEvaluations::value_at(size_t generation) const noexcept
-    {
-        GAPP_ASSERT(generation < data_.size());
-
-        return data_[generation];
-    }
-
 } // namespace gapp::metrics
diff --git a/src/metrics/misc_metrics.hpp b/src/metrics/misc_metrics.hpp
index 7aab1fa5..b640ca93 100644
--- a/src/metrics/misc_metrics.hpp
+++ b/src/metrics/misc_metrics.hpp
@@ -12,9 +12,6 @@ namespace gapp::metrics
     /** Record the number of fitness function evaluations performed in each generation. */
     class FitnessEvaluations final : public Monitor<FitnessEvaluations, std::vector<size_t>>
     {
-    public:
-        size_t value_at(size_t generation) const noexcept;
-    private:
         void initialize(const GaInfo& ga) override;
         void update(const GaInfo& ga) override;
 
diff --git a/src/metrics/monitor.hpp b/src/metrics/monitor.hpp
index 3575977a..d1be4544 100644
--- a/src/metrics/monitor.hpp
+++ b/src/metrics/monitor.hpp
@@ -6,6 +6,8 @@
 #include "monitor_base.hpp"
 #include "../utility/type_id.hpp"
 #include "../utility/type_traits.hpp"
+#include "../utility/concepts.hpp"
+#include "../utility/utility.hpp"
 #include <vector>
 #include <type_traits>
 #include <cstddef>
@@ -17,12 +19,10 @@ namespace gapp::metrics
     * Metrics can be used to track certain attributes of the %GAs
     * in every generation throughout a run.
     * 
-    * New metrics should be derived from this class, and they should implement the following 3 methods:
+    * New metrics should be derived from this class, and they should implement the following methods:
     *   
-    *   - initialize : Used to initialize the monitor at the start of a run.
-    *   - update     : Used to update the monitored attribute once every generation.
-    *   - value_at   : Returns the value of the monitored metric in a specific generation
-                        (must be implemented as a public, non-virtual function).
+    *   - initialize (optional) : Used to initialize the monitor at the start of a run.
+    *   - update                : Used to update the monitored attribute once every generation.
     * 
     * @note The monitor doesn't have access to any encoding specific information about a %GA
     *   by default, so no such information can be tracked.
@@ -30,13 +30,13 @@ namespace gapp::metrics
     * @tparam Derived The type of the derived class.
     * @tparam MetricData The type of the container used to store the monitored metrics (eg. std::vector).
     */
-    template<typename Derived, typename MetricData>
+    template<typename Derived, detail::IndexableContainer MetricData>
     class Monitor : public MonitorBase
     {
     public:
         /** @returns The value of the tracked metric in the specified @p generation. */
         [[nodiscard]]
-        constexpr decltype(auto) operator[](size_t generation) const { return derived().value_at(generation); }
+        constexpr auto operator[](size_t generation) const noexcept { GAPP_ASSERT(generation < data_.size()); return data_[generation]; }
 
         /** @returns The data collected by the monitor throughout the run. */
         [[nodiscard]]
@@ -47,28 +47,17 @@ namespace gapp::metrics
         constexpr size_t size() const noexcept { return data_.size(); }
 
         /** @returns An iterator to the first element of the metric data. */
-        constexpr auto begin() const { return data_.begin(); }
+        constexpr auto begin() const noexcept { return data_.begin(); }
 
         /** @returns An iterator to one past the last element of the metric data. */
-        constexpr auto end() const   { return data_.end(); }
+        constexpr auto end() const noexcept { return data_.end(); }
+
+        void initialize(const GaInfo&) override { data_.clear(); }
 
     protected:
         MetricData data_;
 
-        constexpr Monitor() noexcept
-        {
-            static_assert(!std::is_abstract_v<Derived>,
-                          "The Derived class should implement all the virtual functions of Monitor.");
-            static_assert(detail::is_derived_from_spec_of_v<Derived, Monitor>,
-                          "The first type parameter of the Monitor class must be the derived class.");
-            static_assert(std::is_invocable_v<decltype(&Derived::value_at), const Derived&, size_t>,
-                          "Classes derived from Monitor must implement a public 'value_at(size_t) const' function.");
-        }
-
         size_t type_id() const noexcept final { return detail::type_id<Derived>; }
-
-    private:
-        constexpr const Derived& derived() const noexcept { return static_cast<const Derived&>(*this); }
     };
 
 } // namespace gapp::metrics
diff --git a/src/population/candidate.hpp b/src/population/candidate.hpp
index 5b7fa7f0..5a772f59 100644
--- a/src/population/candidate.hpp
+++ b/src/population/candidate.hpp
@@ -137,27 +137,14 @@ namespace gapp
     template<typename T>
     using CandidatePair = std::pair<Candidate<T>, Candidate<T>>;
 
-    /* Candidates are considered equal if their chromosomes are the same. */
 
+    /**
+    * Comparison operators based on the chromosomes of the candidates.
+    * The comparisons are not transitive if T is a floating-point type. 
+    */
     template<typename T>
     bool operator==(const Candidate<T>& lhs, const Candidate<T>& rhs) noexcept;
 
-    template<typename T>
-    bool operator!=(const Candidate<T>& lhs, const Candidate<T>& rhs) noexcept;
-
-    /* Lexicographical comparison operators based on the chromosomes of the candidates. */
-
-    template<typename T>
-    bool operator<(const Candidate<T>& lhs, const Candidate<T>& rhs) noexcept;
-
-    template<typename T>
-    bool operator<=(const Candidate<T>& lhs, const Candidate<T>& rhs) noexcept;
-
-    template<typename T>
-    bool operator>(const Candidate<T>& lhs, const Candidate<T>& rhs) noexcept;
-
-    template<typename T>
-    bool operator>=(const Candidate<T>& lhs, const Candidate<T>& rhs) noexcept;
 
     /* Hash function for the candidates. */
     template<detail::Hashable T>
@@ -198,45 +185,6 @@ namespace gapp
         }
     }
 
-    template<typename T>
-    bool operator!=(const Candidate<T>& lhs, const Candidate<T>& rhs) noexcept
-    {
-        return !(lhs == rhs);
-    }
-
-    template<typename T>
-    bool operator<(const Candidate<T>& lhs, const Candidate<T>& rhs) noexcept
-    {
-        if constexpr (std::is_floating_point_v<T>)
-        {
-            return std::lexicographical_compare(lhs.chromosome.begin(), lhs.chromosome.end(),
-                                                rhs.chromosome.begin(), rhs.chromosome.end(), 
-                                                math::floatIsLess<T>);
-        }
-        else
-        {
-            return lhs.chromosome < rhs.chromosome;
-        }
-    }
-
-    template<typename T>
-    bool operator>=(const Candidate<T>& lhs, const Candidate<T>& rhs) noexcept
-    {
-        return !(lhs < rhs);
-    }
-
-    template<typename T>
-    bool operator>(const Candidate<T>& lhs, const Candidate<T>& rhs) noexcept
-    {
-        return rhs < lhs;
-    }
-
-    template<typename T>
-    bool operator<=(const Candidate<T>& lhs, const Candidate<T>& rhs) noexcept
-    {
-        return !(rhs < lhs);
-    }
-
     template<detail::Hashable T>
     size_t CandidateHasher<T>::operator()(const Candidate<T>& candidate) const noexcept
     {
diff --git a/src/utility/algorithm.hpp b/src/utility/algorithm.hpp
index dc4d6212..7d5ba2b9 100644
--- a/src/utility/algorithm.hpp
+++ b/src/utility/algorithm.hpp
@@ -43,7 +43,7 @@ namespace gapp::detail
     }
 
     template<std::random_access_iterator Iter, typename Comp = std::less<std::iter_value_t<Iter>>>
-    requires std::strict_weak_order<Comp, std::iter_value_t<Iter>, std::iter_value_t<Iter>>
+    requires std::strict_weak_order<Comp&, std::iter_reference_t<Iter>, std::iter_reference_t<Iter>>
     std::vector<size_t> argsort(Iter first, Iter last, Comp&& comp = {})
     {
         GAPP_ASSERT(std::distance(first, last) >= 0);
@@ -70,7 +70,7 @@ namespace gapp::detail
     }
 
     template<std::random_access_iterator Iter, typename Comp = std::less<std::iter_value_t<Iter>>>
-    requires std::strict_weak_order<Comp, std::iter_value_t<Iter>, std::iter_value_t<Iter>>
+    requires std::strict_weak_order<Comp&, std::iter_reference_t<Iter>, std::iter_reference_t<Iter>>
     std::vector<size_t> partial_argsort(Iter first, Iter middle, Iter last, Comp&& comp = {})
     {
         GAPP_ASSERT(std::distance(first, middle) >= 0);
@@ -104,13 +104,57 @@ namespace gapp::detail
         return indices;
     }
 
-    template<std::random_access_iterator Iter, typename Comp = std::less<std::iter_value_t<Iter>>>
-    requires std::strict_weak_order<Comp, std::iter_value_t<Iter>, std::iter_value_t<Iter>>
-    constexpr size_t argmax(Iter first, Iter last, Comp&& comp = {})
+    template<std::forward_iterator Iter, typename F = std::identity>
+    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);
+        GAPP_ASSERT(std::distance(first, last) >= 0);
+
+        if (first == last) return first;
+
+        Iter max_elem = first;
+        auto max_value = std::invoke(transform, *first);
+
+        for (++first; first != last; ++first)
+        {
+            auto value = std::invoke(transform, *first);
+            if (!(max_value < value)) continue;
+            max_value = std::move(value);
+            max_elem = first;
+        }
+
+        return max_elem;
+    }
+
+    template<std::forward_iterator Iter, typename F = std::identity>
+    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);
+
+        for (++first; first != last; ++first)
+        {
+            auto value = std::invoke(transform, *first);
+            if (!(value < min_value)) continue;
+            min_value = std::move(value);
+            min_elem = first;
+        }
+
+        return min_elem;
+    }
 
-        const auto it = std::max_element(first, last, std::forward<Comp>(comp));
+    template<std::random_access_iterator Iter, typename F = std::identity>
+    requires std::invocable<F&, std::iter_reference_t<Iter>>
+    constexpr size_t argmax(Iter first, Iter last, F&& transform = {})
+    {
+        GAPP_ASSERT(std::distance(first, last) > 0);
+        
+        const auto it = detail::max_element(first, last, std::forward<F>(transform));
         const size_t idx = std::distance(first, it);
 
         if constexpr (detail::is_reverse_iterator_v<Iter>)
@@ -124,13 +168,13 @@ namespace gapp::detail
         }
     }
 
-    template<std::random_access_iterator Iter, typename Comp = std::less<std::iter_value_t<Iter>>>
-    requires std::strict_weak_order<Comp, std::iter_value_t<Iter>, std::iter_value_t<Iter>>
-    constexpr size_t argmin(Iter first, Iter last, Comp&& comp = {})
+    template<std::random_access_iterator Iter, typename F = std::identity>
+    requires std::invocable<F&, std::iter_reference_t<Iter>>
+    constexpr size_t argmin(Iter first, Iter last, F&& transform = {})
     {
         GAPP_ASSERT(std::distance(first, last) > 0);
 
-        const auto it = std::min_element(first, last, std::forward<Comp>(comp));
+        const auto it = detail::min_element(first, last, std::forward<F>(transform));
         const size_t idx = std::distance(first, it);
 
         if constexpr (detail::is_reverse_iterator_v<Iter>)
@@ -145,6 +189,7 @@ namespace gapp::detail
     }
 
     template<std::random_access_iterator Iter, typename URBG>
+    requires std::uniform_random_bit_generator<std::remove_cvref_t<URBG>>
     constexpr void partial_shuffle(Iter first, Iter middle, Iter last, URBG&& gen)
     {
         GAPP_ASSERT(std::distance(first, middle) >= 0);
@@ -160,30 +205,15 @@ namespace gapp::detail
         }
     }
 
-    template<std::input_iterator Iter>
-    constexpr bool contains(Iter first, Iter last, const std::iter_value_t<Iter>& val)
+    template<std::input_iterator Iter, typename T>
+    constexpr bool contains(Iter first, Iter last, const T& val)
     {
         return std::find(first, last, val) != last;
     }
 
-    template<std::forward_iterator Iter, std::predicate<std::iter_value_t<Iter>> Pred>
-    std::vector<Iter> find_all(Iter first, Iter last, Pred&& pred)
-    {
-        GAPP_ASSERT(std::distance(first, last) >= 0);
-
-        std::vector<Iter> result;
-        result.reserve(last - first);
-
-        for (; first != last; ++first)
-        {
-            if (std::invoke(pred, *first)) result.push_back(first);
-        }
-
-        return result;
-    }
-
-    template<std::forward_iterator Iter, std::predicate<std::iter_value_t<Iter>> Pred>
-    auto find_all_v(Iter first, Iter last, Pred&& pred)
+    template<std::forward_iterator Iter, typename Pred>
+    requires std::predicate<Pred&, std::iter_reference_t<Iter>>
+    auto find_all(Iter first, Iter last, Pred&& pred)
     {
         GAPP_ASSERT(std::distance(first, last) >= 0);
 
@@ -200,8 +230,26 @@ namespace gapp::detail
         return result;
     }
 
+    template<detail::IndexableContainer Container, typename T>
+    constexpr std::optional<size_t> index_of(const Container& container, const T& val)
+    {
+        const auto found = std::find(container.begin(), container.end(), val);
+        const size_t idx = std::distance(container.begin(), found);
+
+        return (found != container.end()) ? idx : std::optional<size_t>{};
+    }
 
     template<detail::IndexableContainer Container, std::predicate<typename Container::value_type> Pred>
+    constexpr std::optional<size_t> find_index(const Container& container, Pred&& pred)
+    {
+        const auto found = std::find_if(container.begin(), container.end(), std::forward<Pred>(pred));
+        const size_t idx = std::distance(container.begin(), found);
+
+        return (found != container.end()) ? idx : std::optional<size_t>{};
+    }
+
+    template<detail::IndexableContainer Container, typename Pred>
+    requires std::predicate<Pred&, typename Container::reference>
     std::vector<size_t> find_indices(const Container& container, Pred&& pred)
     {
         std::vector<size_t> indices;
@@ -215,24 +263,6 @@ namespace gapp::detail
         return indices;
     }
 
-    template<detail::IndexableContainer Container>
-    constexpr std::optional<size_t> index_of(const Container& container, const typename Container::value_type& val)
-    {
-        const auto found = std::find(container.begin(), container.end(), val);
-        const size_t idx = std::distance(container.begin(), found);
-        
-        return (found == container.end()) ? std::optional<size_t>{} : idx;
-    }
-
-    template<detail::IndexableContainer Container, std::predicate<typename Container::value_type> Pred>
-    constexpr std::optional<size_t> find_index(const Container& container, Pred&& pred)
-    {
-        const auto found = std::find_if(container.begin(), container.end(), std::forward<Pred>(pred));
-        const size_t idx = std::distance(container.begin(), found);
-
-        return (found == container.end()) ? std::optional<size_t>{} : idx;
-    }
-
     template<typename T>
     std::vector<T> elementwise_min(std::vector<T> left, std::span<const std::type_identity_t<T>> right)
     {
@@ -261,11 +291,11 @@ namespace gapp::detail
         return left;
     }
 
-    template<detail::Container T>
-    constexpr bool erase_first_stable(T& container, const typename T::value_type& value)
+    template<detail::Container Container, typename T>
+    constexpr bool erase_first_stable(Container& container, const T& value)
     {
-        const auto found = std::find(container.cbegin(), container.cend(), value);
-        if (found != container.cend())
+        const auto found = std::find(container.begin(), container.end(), value);
+        if (found != container.end())
         {
             container.erase(found);
             return true;
@@ -277,8 +307,6 @@ namespace gapp::detail
     requires detail::IndexableContainer<std::remove_cvref_t<Container>>
     auto select(Container&& container, const std::vector<size_t>& indices)
     {
-        GAPP_ASSERT(std::all_of(indices.begin(), indices.end(), [&](size_t idx) { return idx < container.size(); }));
-
         using ValueType = typename std::remove_cvref_t<Container>::value_type;
 
         std::vector<ValueType> selected;
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<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))
                 {
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<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();
 
diff --git a/src/utility/functional.hpp b/src/utility/functional.hpp
index e3b6e974..239f4601 100644
--- a/src/utility/functional.hpp
+++ b/src/utility/functional.hpp
@@ -25,61 +25,10 @@ namespace gapp
 
 namespace gapp::detail
 {
-    template<typename T>
-    requires(std::is_lvalue_reference_v<T>)
-    constexpr auto lforward(std::remove_reference_t<T>& t) noexcept
-    {
-        return std::ref<std::remove_reference_t<T>>(t);
-    }
-
-    template<typename T>
-    requires(!std::is_lvalue_reference_v<T>)
-    constexpr T&& lforward(std::remove_reference_t<T>& t) noexcept
-    {
-        return static_cast<T&&>(t);
-    }
-
-    template<typename T>
-    requires(!std::is_lvalue_reference_v<T>)
-    constexpr T&& lforward(std::remove_reference_t<T>&& t) noexcept
-    {
-        return static_cast<T&&>(t);
-    }
-
-    template<typename T>
-    constexpr T drop_rvalue_ref(T&& t) noexcept
-    {
-        return std::forward<T>(t);
-    }
-
-    template<typename F>
-    constexpr auto compose(F&& f) noexcept /* innermost */
-    {
-        return [f = lforward<F>(f)] <typename... Args> (Args&&... args) mutable
-        noexcept(std::is_nothrow_invocable_v<F, Args...>) -> decltype(auto)
-        requires std::is_invocable_v<F, Args...>
-        {
-            return std::invoke(f, std::forward<Args>(args)...);
-        };
-    }
-
-    template<typename F, typename... Fs>
-    constexpr auto compose(F&& f, Fs&&... fs) noexcept
-    {
-        return [f = lforward<F>(f), ...fs = lforward<Fs>(fs)] <typename... Args> (Args&&... args) mutable
-        noexcept(std::is_nothrow_invocable_v<F, Args...> &&
-                 std::is_nothrow_invocable_v<decltype(compose(fs...)), std::invoke_result_t<F, Args...>>) -> decltype(auto)
-        requires(std::is_invocable_v<F, Args...> &&
-                 std::is_invocable_v<decltype(compose(fs...)), std::invoke_result_t<F, Args...>>)
-        {
-            return drop_rvalue_ref(compose(fs...)(std::invoke(f, std::forward<Args>(args)...)));
-        };
-    }
-
     template<typename ValueType, std::invocable<ValueType&> F>
     auto map(const std::vector<ValueType>& cont, F&& f)
     {
-        using MappedType = std::invoke_result_t<std::remove_reference_t<F>, ValueType>;
+        using MappedType = std::invoke_result_t<F&, ValueType>;
         using ResultType = std::vector<std::decay_t<MappedType>>;
 
         ResultType result;
@@ -95,7 +44,7 @@ namespace gapp::detail
     template<typename ValueType, std::invocable<ValueType&> F>
     auto map(const std::vector<ValueType>& cont, F&& f) requires std::is_scalar_v<ValueType>
     {
-        using MappedType = std::invoke_result_t<std::remove_reference_t<F>, ValueType>;
+        using MappedType = std::invoke_result_t<F&, ValueType>;
         using ResultType = std::vector<std::decay_t<MappedType>>;
 
         ResultType result(cont.size());
@@ -108,10 +57,10 @@ namespace gapp::detail
     }
 
     template<typename ValueType, size_t N, typename F>
-    requires std::invocable<std::remove_reference_t<F>, ValueType&>
+    requires std::invocable<F&, ValueType&>
     constexpr auto map(const std::array<ValueType, N>& cont, F&& f)
     {
-        using MappedType = std::invoke_result_t<std::remove_reference_t<F>, ValueType>;
+        using MappedType = std::invoke_result_t<F&, ValueType>;
         using ResultType = std::array<std::decay_t<MappedType>, N>;
 
         ResultType result;
diff --git a/src/utility/math.hpp b/src/utility/math.hpp
index 97035b3f..fa41595e 100644
--- a/src/utility/math.hpp
+++ b/src/utility/math.hpp
@@ -7,6 +7,7 @@
 #include <vector>
 #include <span>
 #include <atomic>
+#include <compare>
 #include <concepts>
 #include <limits>
 #include <cstdint>
@@ -101,33 +102,33 @@ namespace gapp::math
     };
 
 
-    /* Comparison function for floating point numbers. Returns -1 if (lhs < rhs), +1 if (lhs > rhs), and 0 if (lhs == rhs). */
+    /* Three-way comparison function for floating point numbers. */
     template<std::floating_point T>
-    constexpr std::int8_t floatCompare(T lhs, T rhs) noexcept;
+    std::weak_ordering floatCompare(T lhs, T rhs) noexcept;
 
     /* Equality comparison for floating point numbers. Returns true if lhs is approximately equal to rhs. */
     template<std::floating_point T>
-    constexpr bool floatIsEqual(T lhs, T rhs) noexcept;
+    bool floatIsEqual(T lhs, T rhs) noexcept;
 
     /* Less than comparison for floating point numbers. Returns true if lhs is definitely less than rhs. */
     template<std::floating_point T>
-    constexpr bool floatIsLess(T lhs, T rhs) noexcept;
+    bool floatIsLess(T lhs, T rhs) noexcept;
 
     /* Less than comparison for fp numbers. Assumes that lhs is not greater than rhs. */
     template<std::floating_point T>
-    constexpr bool floatIsLessAssumeNotGreater(T lhs, T rhs) noexcept;
+    bool floatIsLessAssumeNotGreater(T lhs, T rhs) noexcept;
 
     /* Greater than comparison for floating point numbers. Returns true if lhs is definitely greater than rhs. */
     template<std::floating_point T>
-    constexpr bool floatIsGreater(T lhs, T rhs) noexcept;
+    bool floatIsGreater(T lhs, T rhs) noexcept;
 
     /* Less than or equal to comparison for floating point numbers. Returns true if lhs is less than or approximately equal to rhs. */
     template<std::floating_point T>
-    constexpr bool floatIsLessEq(T lhs, T rhs) noexcept;
+    bool floatIsLessEq(T lhs, T rhs) noexcept;
 
     /* Greater than or equal to comparison for floating point numbers. Returns true if lhs is greater than or approximately equal to rhs. */
     template<std::floating_point T>
-    constexpr bool floatIsGreaterEq(T lhs, T rhs) noexcept;
+    bool floatIsGreaterEq(T lhs, T rhs) noexcept;
 
     /* Equality comparison for fp vectors. Returns true if the elements of the ranges are approximately equal. */
     template<std::floating_point T>
@@ -182,7 +183,7 @@ namespace gapp::math
 namespace gapp::math
 {
     template<std::floating_point T>
-    constexpr std::int8_t floatCompare(T lhs, T rhs) noexcept
+    std::weak_ordering floatCompare(T lhs, T rhs) noexcept
     {
         GAPP_ASSERT(!std::isnan(lhs) && !std::isnan(rhs));
 
@@ -190,13 +191,13 @@ namespace gapp::math
         const T scale = std::min(std::max(std::abs(lhs), std::abs(rhs)), std::numeric_limits<T>::max());
         const T tol = std::max(Tolerances::rel<T>(scale), Tolerances::abs<T>());
 
-        if (diff >  tol) return  1;  // lhs < rhs
-        if (diff < -tol) return -1;  // lhs > rhs
-        return 0;                    // lhs == rhs
+        if (diff >  tol) return std::weak_ordering::greater;
+        if (diff < -tol) return std::weak_ordering::less;
+        return std::weak_ordering::equivalent;
     }
 
     template<std::floating_point T>
-    constexpr bool floatIsEqual(T lhs, T rhs) noexcept
+    bool floatIsEqual(T lhs, T rhs) noexcept
     {
         const T scale = std::max(std::abs(lhs), std::abs(rhs));
 
@@ -206,7 +207,7 @@ namespace gapp::math
     }
 
     template<std::floating_point T>
-    constexpr bool floatIsLess(T lhs, T rhs) noexcept
+    bool floatIsLess(T lhs, T rhs) noexcept
     {
         const T scale = std::max(std::abs(lhs), std::abs(rhs));
 
@@ -216,7 +217,7 @@ namespace gapp::math
     }
 
     template<std::floating_point T>
-    constexpr bool floatIsLessAssumeNotGreater(T lhs, T rhs) noexcept
+    bool floatIsLessAssumeNotGreater(T lhs, T rhs) noexcept
     {
         const T scale = std::abs(rhs);
 
@@ -226,7 +227,7 @@ namespace gapp::math
     }
 
     template<std::floating_point T>
-    constexpr bool floatIsGreater(T lhs, T rhs) noexcept
+    bool floatIsGreater(T lhs, T rhs) noexcept
     {
         const T scale = std::max(std::abs(lhs), std::abs(rhs));
 
@@ -236,13 +237,13 @@ namespace gapp::math
     }
 
     template<std::floating_point T>
-    constexpr bool floatIsLessEq(T lhs, T rhs) noexcept
+    bool floatIsLessEq(T lhs, T rhs) noexcept
     {
         return !floatIsGreater(lhs, rhs);
     }
 
     template<std::floating_point T>
-    constexpr bool floatIsGreaterEq(T lhs, T rhs) noexcept
+    bool floatIsGreaterEq(T lhs, T rhs) noexcept
     {
         return !floatIsLess(lhs, rhs);
     }
diff --git a/src/utility/type_traits.hpp b/src/utility/type_traits.hpp
index ac1a6cf2..65738683 100644
--- a/src/utility/type_traits.hpp
+++ b/src/utility/type_traits.hpp
@@ -42,31 +42,6 @@ namespace gapp::detail
 
 
 
-    namespace _
-    {
-        template<size_t N, typename... Ts>
-        struct number_of_types_impl : std::integral_constant<size_t, N> {};
-
-        template<size_t N, typename T, typename... Ts>
-        struct number_of_types_impl<N, T, Ts...>
-            : std::integral_constant<size_t, number_of_types_impl<N + 1, Ts...>::value>
-        {};
-
-    } // namespace _
-
-    template<typename... Ts>
-    struct number_of_types : std::integral_constant<size_t, 0> {};
-
-    template<typename T, typename... Ts>
-    struct number_of_types<T, Ts...>
-        : std::integral_constant<size_t, _::number_of_types_impl<1, Ts...>::value>
-    {};
-
-    template<typename... Ts>
-    inline constexpr size_t number_of_types_v = number_of_types<Ts...>::value;
-
-
-
     template<size_t N, typename... Args>
     struct nth_type {};
 
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 <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)); };
+}
diff --git a/test/unit/algorithm.cpp b/test/unit/algorithm.cpp
index 765d5a9c..8648ed54 100644
--- a/test/unit/algorithm.cpp
+++ b/test/unit/algorithm.cpp
@@ -102,6 +102,30 @@ TEST_CASE("partial_argsort", "[algorithm]")
     }
 }
 
+TEST_CASE("max_element", "[algorithm]")
+{
+    const std::vector nums = { 4.0, 0.0, 2.0, 5.0, 1.0 };
+
+    REQUIRE(*detail::max_element(nums.begin(), nums.end()) == 5.0);
+    REQUIRE(*detail::max_element(nums.rbegin(), nums.rend()) == 5.0);
+
+    REQUIRE(detail::max_element(nums.begin(), nums.begin()) == nums.begin());
+
+    REQUIRE(*detail::max_element(nums.begin(), nums.end(), std::negate{}) == 0.0);
+}
+
+TEST_CASE("min_element", "[algorithm]")
+{
+    const std::vector nums = { 4.0, 0.0, 2.0, 5.0, 1.0 };
+
+    REQUIRE(*detail::min_element(nums.begin(), nums.end()) == 0.0);
+    REQUIRE(*detail::min_element(nums.rbegin(), nums.rend()) == 0.0);
+
+    REQUIRE(detail::min_element(nums.begin(), nums.begin()) == nums.begin());
+
+    REQUIRE(*detail::min_element(nums.begin(), nums.end(), std::negate{}) == 5.0);
+}
+
 TEST_CASE("argmax", "[algorithm]")
 {
     const std::vector nums = { 4.0, 0.0, 2.0, 5.0, 1.0 };
@@ -109,7 +133,7 @@ TEST_CASE("argmax", "[algorithm]")
     REQUIRE(detail::argmax(nums.begin(), nums.end()) == 3);
     REQUIRE(detail::argmax(nums.rbegin(), nums.rend()) == 3);
 
-    REQUIRE(detail::argmax(nums.begin(), nums.end(), std::greater<>{}) == 1);
+    REQUIRE(detail::argmax(nums.begin(), nums.end(), std::negate{}) == 1);
 
     REQUIRE(detail::argmax(nums.begin(), nums.begin() + 3) == 0);
     REQUIRE(detail::argmax(nums.begin() + 1, nums.end()) == 2);
@@ -122,7 +146,7 @@ TEST_CASE("argmin", "[algorithm]")
     REQUIRE(detail::argmin(nums.begin(), nums.end()) == 1);
     REQUIRE(detail::argmin(nums.rbegin(), nums.rend()) == 1);
 
-    REQUIRE(detail::argmin(nums.begin(), nums.end(), std::greater<>{}) == 3);
+    REQUIRE(detail::argmin(nums.begin(), nums.end(), std::negate{}) == 3);
 
     REQUIRE(detail::argmin(nums.begin() + 2, nums.end()) == 2);
 }
@@ -166,7 +190,7 @@ TEST_CASE("find_all", "[algorithm]")
     const std::vector nums = { 4, 0, 2, 5, 1 };
 
     const auto odd_nums = detail::find_all(nums.begin(), nums.end(), is_odd);
-    REQUIRE(odd_nums == std::vector{ nums.end() - 2, nums.end() - 1 });
+    REQUIRE(odd_nums == std::vector{ 5, 1 });
 
     const auto big_nums = detail::find_all(nums.begin(), nums.end(), is_big);
     REQUIRE(big_nums.empty());
@@ -175,20 +199,6 @@ TEST_CASE("find_all", "[algorithm]")
     REQUIRE(detail::find_all(nums.begin(), nums.end(), always_false).size() == 0);
 }
 
-TEST_CASE("find_all_v", "[algorithm]")
-{
-    const std::vector nums = { 4, 0, 2, 5, 1 };
-
-    const auto odd_nums = detail::find_all_v(nums.begin(), nums.end(), is_odd);
-    REQUIRE(odd_nums == std::vector{ 5, 1 });
-
-    const auto big_nums = detail::find_all_v(nums.begin(), nums.end(), is_big);
-    REQUIRE(big_nums.empty());
-
-    REQUIRE(detail::find_all_v(nums.begin(), nums.end(), always_true).size()  == nums.size());
-    REQUIRE(detail::find_all_v(nums.begin(), nums.end(), always_false).size() == 0);
-}
-
 TEST_CASE("find_indices", "[algorithm]")
 {
     const std::vector nums = { 4, 0, 2, 5, 1 };
diff --git a/test/unit/functional.cpp b/test/unit/functional.cpp
index 179ec559..9197a646 100644
--- a/test/unit/functional.cpp
+++ b/test/unit/functional.cpp
@@ -16,64 +16,6 @@ constexpr T increment(const T& n) noexcept { return n + 1; }
 using namespace gapp::detail;
 
 
-TEST_CASE("lforward", "[functional]")
-{
-    int n = 2;
-
-    STATIC_REQUIRE(std::is_same_v<int&&, decltype(lforward<int>(n))>);
-    STATIC_REQUIRE(std::is_same_v<int&&, decltype(lforward<int>(1))>);
-    STATIC_REQUIRE(std::is_same_v<int&&, decltype(lforward<int>(std::move(n)))>);
-    STATIC_REQUIRE(std::is_same_v<std::reference_wrapper<int>, decltype(lforward<int&>(n))>);
-}
-
-TEST_CASE("compose", "[functional]")
-{
-    SECTION("functions return by value")
-    {
-        auto f = compose(square<double>, increment<double>, increment<double>);
-        using F = decltype(f);
-
-        STATIC_REQUIRE(std::is_nothrow_invocable_v<F, double>);
-        STATIC_REQUIRE(!std::is_invocable_v<F, void*>);
-        STATIC_REQUIRE(std::is_same_v<std::invoke_result_t<F, double>, double>);
-        STATIC_REQUIRE(std::is_same_v<std::invoke_result_t<F, double&>, double>);
-
-        REQUIRE(f(2.0) == 6.0);
-        REQUIRE(f(5) == 27.0);
-    }
-
-    SECTION("functions return references")
-    {
-        auto f = compose(std::identity{}, std::identity{}, std::identity{});
-        using F = decltype(f);
-
-        STATIC_REQUIRE(std::is_nothrow_invocable_v<F, double>);
-        STATIC_REQUIRE(std::is_invocable_v<F, void*>);
-        STATIC_REQUIRE(std::is_same_v<std::invoke_result_t<F, double>, double>);
-        STATIC_REQUIRE(std::is_same_v<std::invoke_result_t<F, double&>, double&>);
-
-        REQUIRE(f(1.5) == 1.5);
-        REQUIRE(f(nullptr) == nullptr);
-
-        auto in = 1;
-        auto& out = f(in);
-        REQUIRE(&in == &out);
-    }
-
-    SECTION("functions return by value and ref")
-    {
-        auto f = compose(square<double>, std::identity{}, square<double>, std::identity{});
-        using F = decltype(f);
-
-        STATIC_REQUIRE(std::is_nothrow_invocable_v<F, double>);
-        STATIC_REQUIRE(!std::is_invocable_v<F, void*>);
-        STATIC_REQUIRE(std::is_same_v<std::invoke_result_t<F, double>, double>);
-        STATIC_REQUIRE(std::is_same_v<std::invoke_result_t<F, double&>, double>);
-
-        REQUIRE(f(2.0) == 16.0);
-    }
-}
-
 TEST_CASE("map", "[functional]")
 {
     const std::vector nums = { 0.0, 1.2, 5.0, 2.5 };
diff --git a/test/unit/math.cpp b/test/unit/math.cpp
index 28eb8780..e02189de 100644
--- a/test/unit/math.cpp
+++ b/test/unit/math.cpp
@@ -150,38 +150,38 @@ TEST_CASE("fp_compare", "[math]")
         ScopedTolerances _(abs, rel);
         INFO("Relative tolerance eps: " << rel << ", absolute tolerance: " << abs);
 
-        REQUIRE(floatCompare(0.0, 0.0) == 0);
-        REQUIRE(floatCompare(0.0, -0.0) == 0);
-        REQUIRE(floatCompare(-0.0, 0.0) == 0);
-        REQUIRE(floatCompare(-0.0, -0.0) == 0);
+        REQUIRE((floatCompare(0.0, 0.0) == 0));
+        REQUIRE((floatCompare(0.0, -0.0) == 0));
+        REQUIRE((floatCompare(-0.0, 0.0) == 0));
+        REQUIRE((floatCompare(-0.0, -0.0) == 0));
 
-        REQUIRE(floatCompare(4.0, 4.0) == 0);
-        REQUIRE(floatCompare(0.0, 4.0) < 0);
-        REQUIRE(floatCompare(4.0, 0.0) > 0);
+        REQUIRE((floatCompare(4.0, 4.0) == 0));
+        REQUIRE((floatCompare(0.0, 4.0) < 0));
+        REQUIRE((floatCompare(4.0, 0.0) > 0));
 
-        REQUIRE(floatCompare(SMALL, SMALL) == 0);
-        REQUIRE(floatCompare(BIG, BIG) == 0);
-        REQUIRE(floatCompare(INF, INF) == 0);
+        REQUIRE((floatCompare(SMALL, SMALL) == 0));
+        REQUIRE((floatCompare(BIG, BIG) == 0));
+        REQUIRE((floatCompare(INF, INF) == 0));
 
-        REQUIRE(floatCompare(-INF, INF) < 0);
-        REQUIRE(floatCompare(INF, -INF) > 0);
-        REQUIRE(floatCompare(INF, INF) == 0);
-        REQUIRE(floatCompare(-INF, -INF) == 0);
+        REQUIRE((floatCompare(-INF, INF) < 0));
+        REQUIRE((floatCompare(INF, -INF) > 0));
+        REQUIRE((floatCompare(INF, INF) == 0));
+        REQUIRE((floatCompare(-INF, -INF) == 0));
 
-        REQUIRE(floatCompare(0.0, INF) < 0);
-        REQUIRE(floatCompare(INF, 0.0) > 0);
+        REQUIRE((floatCompare(0.0, INF) < 0));
+        REQUIRE((floatCompare(INF, 0.0) > 0));
 
-        REQUIRE(floatCompare(0.0, BIG) < 0);
-        REQUIRE(floatCompare(BIG, 0.0) > 0);
+        REQUIRE((floatCompare(0.0, BIG) < 0));
+        REQUIRE((floatCompare(BIG, 0.0) > 0));
 
-        REQUIRE(floatCompare(SMALL, BIG) < 0);
-        REQUIRE(floatCompare(BIG, SMALL) > 0);
+        REQUIRE((floatCompare(SMALL, BIG) < 0));
+        REQUIRE((floatCompare(BIG, SMALL) > 0));
 
-        REQUIRE(floatCompare(SMALL, INF) < 0);
-        REQUIRE(floatCompare(INF, SMALL) > 0);
+        REQUIRE((floatCompare(SMALL, INF) < 0));
+        REQUIRE((floatCompare(INF, SMALL) > 0));
 
-        REQUIRE(floatCompare(BIG, INF) < 0);
-        REQUIRE(floatCompare(INF, BIG) > 0);
+        REQUIRE((floatCompare(BIG, INF) < 0));
+        REQUIRE((floatCompare(INF, BIG) > 0));
     }
 
     SECTION("is_less_not_greater")
diff --git a/test/unit/metrics.cpp b/test/unit/metrics.cpp
index cf96d2dc..b00e7026 100644
--- a/test/unit/metrics.cpp
+++ b/test/unit/metrics.cpp
@@ -31,9 +31,9 @@ TEMPLATE_TEST_CASE("fitness_metrics", "[metrics]", FitnessMin, FitnessMax, Fitne
     REQUIRE(metric.data().size() == num_gen);
 
     REQUIRE(std::all_of(metric.begin(), metric.end(), detail::is_size(num_obj)));
-    REQUIRE(metric.value_at(4).size() == num_obj);
+    REQUIRE(metric[4].size() == num_obj);
 
-    const auto& val = metric.value_at(7);
+    const auto& val = metric[7];
     REQUIRE(std::all_of(val.begin(), val.end(), detail::equal_to(0.0)));
 }
 
@@ -48,7 +48,7 @@ TEST_CASE("nadir_point_metric", "[metrics]")
     REQUIRE(metric.size() == num_gen);
     REQUIRE(std::all_of(metric.begin(), metric.end(), detail::is_size(num_obj)));
 
-    const auto& val = metric.value_at(5);
+    const auto& val = metric[5];
     REQUIRE(std::all_of(val.begin(), val.end(), detail::equal_to(0.0)));
 }
 
diff --git a/test/unit/type_traits.cpp b/test/unit/type_traits.cpp
index e747f3a3..5efb2771 100644
--- a/test/unit/type_traits.cpp
+++ b/test/unit/type_traits.cpp
@@ -26,14 +26,6 @@ TEST_CASE("is_one_of_templates", "[type_traits]")
     STATIC_REQUIRE(!is_one_of_templates_v<std::vector>);
 }
 
-TEST_CASE("number_of_types", "[type_traits]")
-{
-    STATIC_REQUIRE(number_of_types_v<> == 0);
-    STATIC_REQUIRE(number_of_types_v<void> == 1);
-    STATIC_REQUIRE(number_of_types_v<int, int> == 2);
-    STATIC_REQUIRE(number_of_types_v<void, int, void*, double> == 4);
-}
-
 TEST_CASE("nth_type", "[type_traits]")
 {
     STATIC_REQUIRE(std::is_same_v<nth_type_t<1, int, void>, void>);