diff --git a/.clang-tidy b/.clang-tidy
index 8dbf66c5..d012fa5b 100644
--- a/.clang-tidy
+++ b/.clang-tidy
@@ -29,6 +29,7 @@
   -readability-suspicious-call-argument,
   -readability-uppercase-literal-suffix,
   -readability-braces-around-statements,
+  -readability-qualified-auto,
   -bugprone-easily-swappable-parameters,
   -bugprone-unchecked-optional-access,
   -cppcoreguidelines-special-member-functions,
diff --git a/.github/workflows/analysis.yml b/.github/workflows/analysis.yml
index 600d1a16..7650c618 100644
--- a/.github/workflows/analysis.yml
+++ b/.github/workflows/analysis.yml
@@ -44,7 +44,7 @@ jobs:
             { name: cppcheck,   cmake-flag: CMAKE_CXX_CPPCHECK="cppcheck;--version;--verbose;--report-progress;--enable=all;--error-exitcode=1;--std=c++20;--suppressions-list=../.cppcheck-supressions" }
         ]
         include:
-          - pkgs: clang-18 clang-tools-18 clang-tidy-18 iwyu cppcheck
+          - pkgs: clang-18 clang-tools-18 clang-tidy-18 cppcheck
             cxx: clang++-18
 
     defaults:
diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml
index 76231d0b..6f41a7e6 100644
--- a/.github/workflows/windows.yml
+++ b/.github/workflows/windows.yml
@@ -9,18 +9,13 @@ jobs:
       fail-fast: false
       matrix:
         build-type: [ Release, RelWithDebInfo ]
-        platform: [ x64, Win32 ]
+        platform: [ x64 ]
         generator: [ "Visual Studio 17 2022" ]
         compiler: [ 
             { name: msvc, toolset: v143 }, 
             { name: clang, toolset: ClangCL }
         ]
         build-shared: [ "ON", "OFF" ]
-        exclude:
-          - platform: Win32
-            compiler: ClangCL
-          - platform: Win32
-            build-type: Release
 
     defaults:
       run:
@@ -33,9 +28,6 @@ jobs:
     - name: checkout-repo
       uses: actions/checkout@v4
 
-    #- name: setup-compiler
-    #  run: choco install -y --allow-downgrade ${{ matrix.compiler.pkgs }}
-
     - name: setup-catch
       env:
         CMAKE_GENERATOR: ${{ matrix.generator }}
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 23f095fe..ababc919 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -101,7 +101,7 @@ else() # GNU style compiler interface
     # gcc specific options
     if(CMAKE_CXX_COMPILER_ID MATCHES "GNU")
         # -Warray-bounds, -Wstringop-overflow, -Wstringop-overread are regular false positives since g++-12
-	    set(GAPP_WARN_FLAGS "${GAPP_WARN_FLAGS} -Wlogical-op -Wno-array-bounds -Wno-stringop-overflow -Wno-stringop-overread")
+	    set(GAPP_WARN_FLAGS "${GAPP_WARN_FLAGS} -Wlogical-op -Wno-array-bounds -Wno-stringop-overflow -Wno-stringop-overread -Wno-free-nonheap-object")
     endif()
     # clang specific options
     if(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
diff --git a/src/algorithm/algorithm_base.decl.hpp b/src/algorithm/algorithm_base.decl.hpp
index e89028ba..7bea06a3 100644
--- a/src/algorithm/algorithm_base.decl.hpp
+++ b/src/algorithm/algorithm_base.decl.hpp
@@ -1,12 +1,12 @@
 /* Copyright (c) 2022 Krisztián Rugási. Subject to the MIT License. */
 
-#ifndef GA_ALGORITHM_ALGORITHM_BASE_DECL_HPP
-#define GA_ALGORITHM_ALGORITHM_BASE_DECL_HPP
+#ifndef GAPP_ALGORITHM_ALGORITHM_BASE_DECL_HPP
+#define GAPP_ALGORITHM_ALGORITHM_BASE_DECL_HPP
 
 #include "selection_base.hpp"
 #include "replacement_base.hpp"
 #include "../core/population.hpp"
-#include <vector>
+#include "../utility/small_vector.hpp"
 #include <cstddef>
 
 namespace gapp
@@ -138,9 +138,9 @@ namespace gapp::algorithm
         * @param ga The %GA that uses the algorithm.
         * @returns The indices of the pareto optimal solutions in the current population.
         */
-        virtual std::vector<size_t> optimalSolutionsImpl(const GaInfo&) const { return {}; }
+        virtual small_vector<size_t> optimalSolutionsImpl(const GaInfo&) const { return {}; }
     };
 
 } // namespace gapp::algorithm
 
-#endif // !GA_ALGORITHM_ALGORITHM_BASE_DECL_HPP
\ No newline at end of file
+#endif // !GAPP_ALGORITHM_ALGORITHM_BASE_DECL_HPP
\ No newline at end of file
diff --git a/src/algorithm/algorithm_base.hpp b/src/algorithm/algorithm_base.hpp
index f455bc83..085d086f 100644
--- a/src/algorithm/algorithm_base.hpp
+++ b/src/algorithm/algorithm_base.hpp
@@ -1,9 +1,9 @@
 /* Copyright (c) 2022 Krisztián Rugási. Subject to the MIT License. */
 
-#ifndef GA_ALGORITHM_ALGORITHM_BASE_HPP
-#define GA_ALGORITHM_ALGORITHM_BASE_HPP
+#ifndef GAPP_ALGORITHM_ALGORITHM_BASE_HPP
+#define GAPP_ALGORITHM_ALGORITHM_BASE_HPP
 
 #include "algorithm_base.decl.hpp"
 #include "algorithm_base.impl.hpp"
 
-#endif //!GA_ALGORITHM_ALGORITHM_BASE_HPP
\ No newline at end of file
+#endif //!GAPP_ALGORITHM_ALGORITHM_BASE_HPP
\ No newline at end of file
diff --git a/src/algorithm/algorithm_base.impl.hpp b/src/algorithm/algorithm_base.impl.hpp
index fc475f4d..4e15eda1 100644
--- a/src/algorithm/algorithm_base.impl.hpp
+++ b/src/algorithm/algorithm_base.impl.hpp
@@ -1,7 +1,7 @@
 /* Copyright (c) 2022 Krisztián Rugási. Subject to the MIT License. */
 
-#ifndef GA_ALGORITHM_ALGORITHM_BASE_IMPL_HPP
-#define GA_ALGORITHM_ALGORITHM_BASE_IMPL_HPP
+#ifndef GAPP_ALGORITHM_ALGORITHM_BASE_IMPL_HPP
+#define GAPP_ALGORITHM_ALGORITHM_BASE_IMPL_HPP
 
 #include "algorithm_base.decl.hpp"
 #include "../core/ga_info.hpp"
@@ -17,7 +17,7 @@
 namespace gapp::algorithm
 {
     template<typename T>
-    auto Algorithm::select(const GA<T>& ga, const Population<T>& pop, const FitnessMatrix& fmat) const -> const Candidate<T>&
+    const Candidate<T>& Algorithm::select(const GA<T>& ga, const Population<T>& pop, const FitnessMatrix& fmat) const
     {
         GAPP_ASSERT(ga.population_size() == pop.size());
         GAPP_ASSERT(pop.size() == fmat.size());
@@ -30,7 +30,7 @@ namespace gapp::algorithm
     }
 
     template<typename T>
-    auto Algorithm::nextPopulation(const GA<T>& ga, Population<T> parents, Population<T> children) -> Population<T>
+    Population<T> Algorithm::nextPopulation(const GA<T>& ga, Population<T> parents, Population<T> children)
     {
         GAPP_ASSERT(ga.population_size() == parents.size());
         GAPP_ASSERT(ga.population_size() <= children.size());
@@ -48,7 +48,7 @@ namespace gapp::algorithm
     }
 
     template<typename T>
-    auto Algorithm::optimalSolutions(const GA<T>& ga, const Population<T>& pop) const -> Candidates<T>
+    Candidates<T> Algorithm::optimalSolutions(const GA<T>& ga, const Population<T>& pop) const
     {
         GAPP_ASSERT(ga.population_size() == pop.size());
 
@@ -66,4 +66,4 @@ namespace gapp::algorithm
 
 } // namespace gapp::algorithm
 
-#endif //!GA_ALGORITHM_ALGORITHM_BASE_IMPL_HPP
+#endif //!GAPP_ALGORITHM_ALGORITHM_BASE_IMPL_HPP
diff --git a/src/algorithm/nd_sort.cpp b/src/algorithm/nd_sort.cpp
index 37df5156..3f411927 100644
--- a/src/algorithm/nd_sort.cpp
+++ b/src/algorithm/nd_sort.cpp
@@ -278,7 +278,7 @@ namespace gapp::algorithm::dtl
         pareto_fronts.reserve(popsize);
 
         size_t current_rank = 0;
-        std::vector<size_t> removed_rows;
+        small_vector<size_t, 8> removed_rows;
 
         while (pareto_fronts.size() != popsize)
         {
diff --git a/src/algorithm/nsga2.cpp b/src/algorithm/nsga2.cpp
index 1be29687..0843784b 100644
--- a/src/algorithm/nsga2.cpp
+++ b/src/algorithm/nsga2.cpp
@@ -26,11 +26,11 @@ namespace gapp::algorithm
         GAPP_ASSERT(!fmat.empty());
         GAPP_ASSERT(std::none_of(pareto_fronts.begin(), pareto_fronts.end(), detail::is_size(0)));
 
-        std::vector<double> crowding_distances(fmat.size(), 0.0);
+        std::vector crowding_distances(fmat.size(), 0.0);
 
         for (size_t obj = 0; obj < fmat.ncols(); obj++)
         {
-            const FitnessVector fvec = fmat.column(obj);
+            const auto fvec = fmat.column(obj);
 
             for (const auto& front : pareto_fronts)
             {
@@ -83,7 +83,7 @@ namespace gapp::algorithm
         return idx2;
     }
 
-    std::vector<size_t> NSGA2::nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat)
+    small_vector<size_t> NSGA2::nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat)
     {
         GAPP_ASSERT(ga.num_objectives() > 1);
         GAPP_ASSERT(fmat.ncols() == ga.num_objectives());
@@ -103,7 +103,7 @@ namespace gapp::algorithm
         dists_ = crowdingDistances(fmat, pareto_fronts.fronts());
         dists_.resize(popsize);
 
-        std::vector<size_t> new_pop(popsize);
+        small_vector<size_t> new_pop(popsize);
         for (size_t i = 0; i < popsize; i++)
         {
             new_pop[i] = pareto_fronts[i].idx;
@@ -113,7 +113,7 @@ namespace gapp::algorithm
         return new_pop;
     }
 
-    std::vector<size_t> NSGA2::optimalSolutionsImpl(const GaInfo&) const
+    small_vector<size_t> NSGA2::optimalSolutionsImpl(const GaInfo&) const
     {
         return detail::find_indices(ranks_, detail::equal_to(0_sz));
     }
diff --git a/src/algorithm/nsga2.hpp b/src/algorithm/nsga2.hpp
index e1e29219..3bd1d6e8 100644
--- a/src/algorithm/nsga2.hpp
+++ b/src/algorithm/nsga2.hpp
@@ -1,9 +1,10 @@
 /* Copyright (c) 2022 Krisztián Rugási. Subject to the MIT License. */
 
-#ifndef GA_ALGORITHM_NSGA2_HPP
-#define GA_ALGORITHM_NSGA2_HPP
+#ifndef GAPP_ALGORITHM_NSGA2_HPP
+#define GAPP_ALGORITHM_NSGA2_HPP
 
 #include "algorithm_base.hpp"
+#include "../utility/small_vector.hpp"
 #include <vector>
 #include <cstddef>
 
@@ -38,9 +39,9 @@ namespace gapp::algorithm
         void prepareSelectionsImpl(const GaInfo&, const FitnessMatrix&) override {}
         size_t selectImpl(const GaInfo& ga, const FitnessMatrix& fmat) const override;
 
-        std::vector<size_t> nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat) override;
+        small_vector<size_t> nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat) override;
 
-        std::vector<size_t> optimalSolutionsImpl(const GaInfo& ga) const override;
+        small_vector<size_t> optimalSolutionsImpl(const GaInfo& ga) const override;
 
         std::vector<size_t> ranks_;
         std::vector<double> dists_;
@@ -48,4 +49,4 @@ namespace gapp::algorithm
 
 } // namespace gapp::algorithm
 
-#endif // !GA_ALGORITHM_NSGA2_HPP
\ No newline at end of file
+#endif // !GAPP_ALGORITHM_NSGA2_HPP
\ No newline at end of file
diff --git a/src/algorithm/nsga3.cpp b/src/algorithm/nsga3.cpp
index 481f6389..e40f7217 100644
--- a/src/algorithm/nsga3.cpp
+++ b/src/algorithm/nsga3.cpp
@@ -8,6 +8,7 @@
 #include "../metrics/pop_stats.hpp"
 #include "../utility/algorithm.hpp"
 #include "../utility/functional.hpp"
+#include "../utility/small_vector.hpp"
 #include "../utility/thread_pool.hpp"
 #include "../utility/math.hpp"
 #include "../utility/rng.hpp"
@@ -42,11 +43,11 @@ namespace gapp::algorithm
     }
 
     /* Create a weight vector for the given axis (used in the ASF). */
-    static inline std::vector<double> weightVector(size_t dimensions, size_t axis)
+    static inline small_vector<double> weightVector(size_t dimensions, size_t axis)
     {
         GAPP_ASSERT(dimensions > axis);
 
-        std::vector weights(dimensions, 1E-6);
+        small_vector weights(dimensions, 1E-6);
         weights[axis] = 1.0;
 
         return weights;
@@ -127,7 +128,7 @@ namespace gapp::algorithm
         void incrementNicheCount(std::vector<size_t>& refs, size_t ref);
 
         /* Create a new population from pareto_fronts. */
-        std::vector<size_t> createPopulation(std::span<const FrontElement> pareto_fronts);
+        small_vector<size_t> createPopulation(std::span<const FrontElement> pareto_fronts);
     };
 
 
@@ -296,9 +297,9 @@ namespace gapp::algorithm
         std::iter_swap(current, std::prev(first_eq));
     }
 
-    std::vector<size_t> NSGA3::Impl::createPopulation(std::span<const FrontElement> pareto_fronts)
+    small_vector<size_t> NSGA3::Impl::createPopulation(std::span<const FrontElement> pareto_fronts)
     {
-        std::vector<size_t> new_pop;
+        small_vector<size_t> new_pop;
         std::vector<Impl::CandidateInfo> new_info;
 
         new_pop.reserve(pareto_fronts.size());
@@ -340,7 +341,7 @@ namespace gapp::algorithm
         pimpl_->recalcNicheCounts(pareto_fronts);
     }
 
-    std::vector<size_t> NSGA3::nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat)
+    small_vector<size_t> NSGA3::nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat)
     {
         GAPP_ASSERT(ga.num_objectives() > 1);
         GAPP_ASSERT(fmat.ncols() == ga.num_objectives());
@@ -403,7 +404,7 @@ namespace gapp::algorithm
         return pimpl_->nichedCompare(idx1, idx2) ? idx1 : idx2;
     }
 
-    std::vector<size_t> NSGA3::optimalSolutionsImpl(const GaInfo&) const
+    small_vector<size_t> NSGA3::optimalSolutionsImpl(const GaInfo&) const
     {
         return detail::find_indices(pimpl_->sol_info_, [](const Impl::CandidateInfo& sol) { return sol.rank == 0; });
     }
diff --git a/src/algorithm/nsga3.hpp b/src/algorithm/nsga3.hpp
index 24aaea07..f115b8fb 100644
--- a/src/algorithm/nsga3.hpp
+++ b/src/algorithm/nsga3.hpp
@@ -1,12 +1,12 @@
 /* Copyright (c) 2022 Krisztián Rugási. Subject to the MIT License. */
 
-#ifndef GA_ALGORITHM_NSGA3_HPP
-#define GA_ALGORITHM_NSGA3_HPP
+#ifndef GAPP_ALGORITHM_NSGA3_HPP
+#define GAPP_ALGORITHM_NSGA3_HPP
 
 #include "algorithm_base.hpp"
 #include "reference_lines.hpp"
 #include "../utility/math.hpp"
-#include <vector>
+#include "../utility/small_vector.hpp"
 #include <functional>
 #include <memory>
 #include <cstddef>
@@ -64,9 +64,9 @@ namespace gapp::algorithm
         void initializeImpl(const GaInfo& ga) override;
         size_t selectImpl(const GaInfo& ga, const FitnessMatrix& fmat) const override;
 
-        std::vector<size_t> nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat) override;
+        small_vector<size_t> nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat) override;
 
-        std::vector<size_t> optimalSolutionsImpl(const GaInfo& ga) const override;
+        small_vector<size_t> optimalSolutionsImpl(const GaInfo& ga) const override;
 
         struct Impl;
         std::unique_ptr<Impl> pimpl_;
@@ -74,4 +74,4 @@ namespace gapp::algorithm
 
 } // namespace gapp::algorithm
 
-#endif // !GA_ALGORITHM_NSGA2_HPP
\ No newline at end of file
+#endif // !GAPP_ALGORITHM_NSGA2_HPP
\ No newline at end of file
diff --git a/src/algorithm/replacement_base.hpp b/src/algorithm/replacement_base.hpp
index 56a64132..4ef23ed1 100644
--- a/src/algorithm/replacement_base.hpp
+++ b/src/algorithm/replacement_base.hpp
@@ -1,10 +1,10 @@
 /* Copyright (c) 2023 Krisztián Rugási. Subject to the MIT License. */
 
-#ifndef GA_ALGORITHM_SOGA_REPLACEMENT_BASE_HPP
-#define GA_ALGORITHM_SOGA_REPLACEMENT_BASE_HPP
+#ifndef GAPP_ALGORITHM_SOGA_REPLACEMENT_BASE_HPP
+#define GAPP_ALGORITHM_SOGA_REPLACEMENT_BASE_HPP
 
 #include "../core/population.hpp"
-#include <vector>
+#include "../utility/small_vector.hpp"
 #include <cstddef>
 
 namespace gapp
@@ -44,7 +44,7 @@ namespace gapp::replacement
         * @param fmat The fitness matrix of the combined parent and child populations.
         * @returns The indices of the candidates selected from the fitness matrix.
         */
-        virtual std::vector<size_t> nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat) = 0;
+        virtual small_vector<size_t> nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat) = 0;
         
 
         /** Destructor. */
@@ -61,4 +61,4 @@ namespace gapp::replacement
 
 } // namespace gapp::replacement
 
-#endif // !GA_ALGORITHM_SOGA_REPLACEMENT_BASE_HPP
\ No newline at end of file
+#endif // !GAPP_ALGORITHM_SOGA_REPLACEMENT_BASE_HPP
\ No newline at end of file
diff --git a/src/algorithm/single_objective.cpp b/src/algorithm/single_objective.cpp
index 8c3a5f63..491b0cc2 100644
--- a/src/algorithm/single_objective.cpp
+++ b/src/algorithm/single_objective.cpp
@@ -59,7 +59,7 @@ namespace gapp::algorithm
         selection_->initializeImpl(ga);
     }
 
-    std::vector<size_t> SingleObjective::nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat)
+    small_vector<size_t> SingleObjective::nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat)
     {
         GAPP_ASSERT(replacement_);
         GAPP_ASSERT(ga.num_objectives() == 1, "The number of objectives must be 1 for the single-objective algorithms.");
diff --git a/src/algorithm/single_objective.hpp b/src/algorithm/single_objective.hpp
index 119ae469..a7e99595 100644
--- a/src/algorithm/single_objective.hpp
+++ b/src/algorithm/single_objective.hpp
@@ -1,15 +1,15 @@
 /* Copyright (c) 2022 Krisztián Rugási. Subject to the MIT License. */
 
-#ifndef GA_ALGORITHM_SINGLE_OBJECTIVE_HPP
-#define GA_ALGORITHM_SINGLE_OBJECTIVE_HPP
+#ifndef GAPP_ALGORITHM_SINGLE_OBJECTIVE_HPP
+#define GAPP_ALGORITHM_SINGLE_OBJECTIVE_HPP
 
 #include "algorithm_base.decl.hpp"
 #include "selection_base.hpp"
 #include "replacement_base.hpp"
 #include "soga_selection.hpp"
 #include "soga_replacement.hpp"
+#include "../utility/small_vector.hpp"
 #include "../utility/utility.hpp"
-#include <vector>
 #include <concepts>
 #include <functional>
 #include <memory>
@@ -45,7 +45,7 @@ namespace gapp::algorithm
         * when not using a replacement policy derived from replacement::Replacement.
         * @see replacement_method()
         */
-        using ReplacementCallable = std::function<std::vector<size_t>(const GaInfo&, const FitnessMatrix&)>;
+        using ReplacementCallable = std::function<small_vector<size_t>(const GaInfo&, const FitnessMatrix&)>;
 
 
         /**
@@ -153,7 +153,7 @@ namespace gapp::algorithm
         void prepareSelectionsImpl(const GaInfo& ga, const FitnessMatrix& fmat) override;
         size_t selectImpl(const GaInfo& ga, const FitnessMatrix& fmat) const override;
 
-        std::vector<size_t> nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat) override;
+        small_vector<size_t> nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat) override;
 
         std::unique_ptr<selection::Selection> selection_;
         std::unique_ptr<replacement::Replacement> replacement_;
@@ -204,4 +204,4 @@ namespace gapp::algorithm
 
 } // namespace gapp::algorithm
 
-#endif // !GA_ALGORITHM_SINGLE_OBJECTIVE_HPP
\ No newline at end of file
+#endif // !GAPP_ALGORITHM_SINGLE_OBJECTIVE_HPP
\ No newline at end of file
diff --git a/src/algorithm/soga_replacement.cpp b/src/algorithm/soga_replacement.cpp
index 92d24483..1afe23b1 100644
--- a/src/algorithm/soga_replacement.cpp
+++ b/src/algorithm/soga_replacement.cpp
@@ -12,12 +12,12 @@
 
 namespace gapp::replacement
 {
-    std::vector<size_t> KeepChildren::nextPopulationImpl(const GaInfo& ga, const FitnessMatrix&)
+    small_vector<size_t> KeepChildren::nextPopulationImpl(const GaInfo& ga, const FitnessMatrix&)
     {
         return detail::index_vector(ga.population_size(), ga.population_size());
     }
 
-    std::vector<size_t> Elitism::nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat)
+    small_vector<size_t> Elitism::nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat)
     {
         GAPP_ASSERT(fmat.size() >= 2 * ga.population_size());
 
@@ -29,14 +29,14 @@ namespace gapp::replacement
             return math::paretoCompareLess(rhs, lhs); // descending
         });
 
-        std::vector<size_t> indices(ga.population_size());
+        small_vector<size_t> indices(ga.population_size());
         std::copy(sorted_parent_indices.begin(), sorted_parent_indices.begin() + elite_count, indices.begin());
         std::iota(indices.begin() + elite_count, indices.end(), ga.population_size());
 
         return indices;
     }
 
-    std::vector<size_t> KeepBest::nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat)
+    small_vector<size_t> KeepBest::nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat)
     {
         GAPP_ASSERT(fmat.size() >= ga.population_size());
 
@@ -58,7 +58,7 @@ namespace gapp::replacement
         replacement_ = std::move(f);
     }
 
-    std::vector<size_t> Lambda::nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat)
+    small_vector<size_t> Lambda::nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat)
     {
         GAPP_ASSERT(replacement_);
 
diff --git a/src/algorithm/soga_replacement.hpp b/src/algorithm/soga_replacement.hpp
index bed29b59..974490ee 100644
--- a/src/algorithm/soga_replacement.hpp
+++ b/src/algorithm/soga_replacement.hpp
@@ -1,11 +1,11 @@
 /* Copyright (c) 2022 Krisztián Rugási. Subject to the MIT License. */
 
-#ifndef GA_ALGORITHM_SOGA_REPLACEMENT_HPP
-#define GA_ALGORITHM_SOGA_REPLACEMENT_HPP
+#ifndef GAPP_ALGORITHM_SOGA_REPLACEMENT_HPP
+#define GAPP_ALGORITHM_SOGA_REPLACEMENT_HPP
 
 #include "replacement_base.hpp"
 #include "../core/population.hpp"
-#include <vector>
+#include "../utility/small_vector.hpp"
 #include <functional>
 #include <cstddef>
 
@@ -28,7 +28,7 @@ namespace gapp::replacement
     class KeepChildren final : public Replacement
     {
     private:
-        std::vector<size_t> nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat) override;
+        small_vector<size_t> nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat) override;
     };
 
 
@@ -69,7 +69,7 @@ namespace gapp::replacement
         constexpr size_t elite_num() const noexcept { return n_; }
 
     private:
-        std::vector<size_t> nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat) override;
+        small_vector<size_t> nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat) override;
 
         size_t n_;
     };
@@ -85,7 +85,7 @@ namespace gapp::replacement
     class KeepBest final : public Replacement
     {
     private:
-        std::vector<size_t> nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat) override;
+        small_vector<size_t> nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat) override;
     };
 
 
@@ -96,12 +96,12 @@ namespace gapp::replacement
     class Lambda final : public Replacement
     {
     public:
-        using ReplacementCallable = std::function<std::vector<size_t>(const GaInfo&, const FitnessMatrix&)>;
+        using ReplacementCallable = std::function<small_vector<size_t>(const GaInfo&, const FitnessMatrix&)>;
 
         explicit Lambda(ReplacementCallable f) noexcept;
 
     private:
-        std::vector<size_t> nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat) override;
+        small_vector<size_t> nextPopulationImpl(const GaInfo& ga, const FitnessMatrix& fmat) override;
 
         ReplacementCallable replacement_;
     };
diff --git a/src/core/candidate.hpp b/src/core/candidate.hpp
index 69476b6b..7ba1d5d3 100644
--- a/src/core/candidate.hpp
+++ b/src/core/candidate.hpp
@@ -1,9 +1,10 @@
 /* Copyright (c) 2022 Krisztián Rugási. Subject to the MIT License. */
 
-#ifndef GA_CORE_CANDIDATE_HPP
-#define GA_CORE_CANDIDATE_HPP
+#ifndef GAPP_CORE_CANDIDATE_HPP
+#define GAPP_CORE_CANDIDATE_HPP
 
 #include "../utility/math.hpp"
+#include "../utility/small_vector.hpp"
 #include "../utility/matrix.hpp"
 #include "../utility/concepts.hpp"
 #include <vector>
@@ -14,7 +15,7 @@
 namespace gapp
 {
     /** The class used to represent the fitness of the candidates. Contains a fitness value for each objective. */
-    using FitnessVector = std::vector<double>;
+    using FitnessVector = small_vector<double>;
 
     /** 
     * The class used to represent the fitness values of multiple candidates.
@@ -128,8 +129,8 @@ namespace gapp
         Candidate& operator=(Candidate&&)       = default;
         ~Candidate()                            = default;
 
-        Chromosome<T> chromosome;   /**< The chromosome encoding the solution. */
         FitnessVector fitness;      /**< The fitness values of the solution (for every objective). */
+        Chromosome<T> chromosome;   /**< The chromosome encoding the solution. */
         bool is_evaluated = false;  /**< True if the candidate's fitness value doesn't need to be computed. */
     };
 
@@ -215,4 +216,4 @@ namespace std
 
 } // namespace std
 
-#endif // !GA_CORE_CANDIDATE_HPP
\ No newline at end of file
+#endif // !GAPP_CORE_CANDIDATE_HPP
\ No newline at end of file
diff --git a/src/core/population.cpp b/src/core/population.cpp
index c7a716a8..9f84a92f 100644
--- a/src/core/population.cpp
+++ b/src/core/population.cpp
@@ -3,6 +3,7 @@
 #include "population.hpp"
 #include "../utility/math.hpp"
 #include "../utility/algorithm.hpp"
+#include "../utility/functional.hpp"
 #include <algorithm>
 #include <iterator>
 #include <vector>
@@ -13,19 +14,19 @@ namespace gapp::detail
     FitnessVector toFitnessVector(FitnessMatrix::const_iterator first, FitnessMatrix::const_iterator last)
     {
         FitnessVector fitness_vector(last - first);
-        std::transform(first, last, fitness_vector.begin(), [](const auto& row) { return row[0]; });
+        std::transform(first, last, fitness_vector.begin(), detail::element_at(0));
 
         return fitness_vector;
     }
 
-    std::vector<size_t> findParetoFront(const FitnessMatrix& fmat)
+    small_vector<size_t> findParetoFront(const FitnessMatrix& fmat)
     {
         if (fmat.empty()) return {};
 
         return (fmat.ncols() == 1) ? findParetoFront1D(fmat) : findParetoFrontBest(fmat);
     }
 
-    std::vector<size_t> findParetoFront1D(const FitnessMatrix& fmat)
+    small_vector<size_t> findParetoFront1D(const FitnessMatrix& fmat)
     {
         const auto best = std::max_element(fmat.begin(), fmat.end(),
         [](const auto& lhs, const auto& rhs) noexcept
@@ -41,7 +42,7 @@ namespace gapp::detail
         });
     }
 
-    std::vector<size_t> findParetoFrontSort(const FitnessMatrix& fmat)
+    small_vector<size_t> findParetoFrontSort(const FitnessMatrix& fmat)
     {
         const auto indices = detail::argsort(fmat.begin(), fmat.end(), [](const auto& lhs, const auto& rhs) noexcept
         {
@@ -52,7 +53,7 @@ namespace gapp::detail
             return false;
         });
 
-        std::vector<size_t> optimal_indices;
+        small_vector<size_t> optimal_indices;
 
         for (size_t idx : indices)
         {
@@ -67,14 +68,14 @@ namespace gapp::detail
         return optimal_indices;
     }
 
-    std::vector<size_t> findParetoFrontBest(const FitnessMatrix& fmat)
+    small_vector<size_t> findParetoFrontBest(const FitnessMatrix& fmat)
     {
         // Implementation of the BEST algorithm based on the description in:
         //      Godfrey et al. "Algorithms and analyses for maximal vector computation." The VLDB Journal 16, no. 1 (2007): 5-28.
 
         auto indices = detail::index_vector(fmat.size());
 
-        std::vector<size_t> optimal_indices;
+        small_vector<size_t> optimal_indices;
         optimal_indices.reserve(fmat.size());
 
         auto first = indices.begin();
@@ -128,7 +129,7 @@ namespace gapp::detail
         return optimal_indices;
     }
 
-    std::vector<size_t> findParetoFrontKungImpl(const FitnessMatrix& fmat, std::vector<size_t>::const_iterator first, std::vector<size_t>::const_iterator last)
+    small_vector<size_t> findParetoFrontKungImpl(const FitnessMatrix& fmat, small_vector<size_t>::const_iterator first, small_vector<size_t>::const_iterator last)
     {
         if (std::distance(first, last) == 1) return { *first };
 
@@ -154,7 +155,7 @@ namespace gapp::detail
         return top_half;
     }
 
-    std::vector<size_t> findParetoFrontKung(const FitnessMatrix& fmat)
+    small_vector<size_t> findParetoFrontKung(const FitnessMatrix& fmat)
     {
         /* See: Kung et al. "On finding the maxima of a set of vectors." Journal of the ACM (JACM) 22.4 (1975): 469-476. */
         /* Doesn't work for d = 1 (single-objective optimization). */
diff --git a/src/core/population.hpp b/src/core/population.hpp
index e3844c63..46436b57 100644
--- a/src/core/population.hpp
+++ b/src/core/population.hpp
@@ -1,10 +1,10 @@
 /* Copyright (c) 2022 Krisztián Rugási. Subject to the MIT License. */
 
-#ifndef GA_CORE_POPULATION_HPP
-#define GA_CORE_POPULATION_HPP
+#ifndef GAPP_CORE_POPULATION_HPP
+#define GAPP_CORE_POPULATION_HPP
 
 #include "candidate.hpp"
-#include <vector>
+#include "../utility/small_vector.hpp"
 #include <cstddef>
 
 namespace gapp
@@ -32,12 +32,12 @@ namespace gapp::detail
     template<typename T>
     Candidates<T> findParetoFront(const Population<T>& pop);
 
-    std::vector<size_t> findParetoFront(const FitnessMatrix& fmat);
+    small_vector<size_t> findParetoFront(const FitnessMatrix& fmat);
 
-    std::vector<size_t> findParetoFront1D(const FitnessMatrix& fmat);
-    std::vector<size_t> findParetoFrontSort(const FitnessMatrix& fmat);
-    std::vector<size_t> findParetoFrontBest(const FitnessMatrix& fmat);
-    std::vector<size_t> findParetoFrontKung(const FitnessMatrix& fmat);
+    small_vector<size_t> findParetoFront1D(const FitnessMatrix& fmat);
+    small_vector<size_t> findParetoFrontSort(const FitnessMatrix& fmat);
+    small_vector<size_t> findParetoFrontBest(const FitnessMatrix& fmat);
+    small_vector<size_t> findParetoFrontKung(const FitnessMatrix& fmat);
 
     /* Find the pareto-optimal solutions in the set (lhs U rhs), assuming both lhs and rhs are pareto sets. */
     template<typename T>
@@ -169,4 +169,4 @@ namespace gapp::detail
 
 } // namespace gapp::detail
 
-#endif // !GA_CORE_POPULATION_HPP
\ No newline at end of file
+#endif // !GAPP_CORE_POPULATION_HPP
\ No newline at end of file
diff --git a/src/crossover/neighbour_list.hpp b/src/crossover/neighbour_list.hpp
index 7e492436..37298090 100644
--- a/src/crossover/neighbour_list.hpp
+++ b/src/crossover/neighbour_list.hpp
@@ -27,10 +27,10 @@ namespace gapp::crossover::dtl
         using size_type              = std::size_t;
         using difference_type        = std::ptrdiff_t;
 
-        using iterator               = typename std::vector<value_type>::iterator;
-        using const_iterator         = typename std::vector<value_type>::const_iterator;
-        using reverse_iterator       = typename std::vector<value_type>::reverse_iterator;
-        using const_reverse_iterator = typename std::vector<value_type>::const_reverse_iterator;
+        using iterator               = typename small_vector<value_type>::iterator;
+        using const_iterator         = typename small_vector<value_type>::const_iterator;
+        using reverse_iterator       = typename small_vector<value_type>::reverse_iterator;
+        using const_reverse_iterator = typename small_vector<value_type>::const_reverse_iterator;
 
         constexpr auto begin() noexcept { return neighbours_.begin(); }
         constexpr auto end() noexcept { return neighbours_.end(); }
diff --git a/src/metrics/metric_set.hpp b/src/metrics/metric_set.hpp
index 6d34bb55..bd4962d4 100644
--- a/src/metrics/metric_set.hpp
+++ b/src/metrics/metric_set.hpp
@@ -4,7 +4,7 @@
 #define GA_METRICS_METRIC_SET_HPP
 
 #include "monitor_base.hpp"
-#include <vector>
+#include "../utility/small_vector.hpp"
 #include <memory>
 #include <type_traits>
 #include <concepts>
@@ -33,7 +33,7 @@ namespace gapp::detail
         void update(const GaInfo& ga);
 
     private:
-        std::vector<std::unique_ptr<metrics::MonitorBase>> metrics_;
+        small_vector<std::unique_ptr<metrics::MonitorBase>> metrics_;
     };
 
 } // namespace gapp::detail
diff --git a/src/utility/algorithm.hpp b/src/utility/algorithm.hpp
index bee34fae..5ca02af7 100644
--- a/src/utility/algorithm.hpp
+++ b/src/utility/algorithm.hpp
@@ -1,8 +1,9 @@
 /* Copyright (c) 2022 Krisztián Rugási. Subject to the MIT License. */
 
-#ifndef GA_UTILITY_ALGORITHM_HPP
-#define GA_UTILITY_ALGORITHM_HPP
+#ifndef GAPP_UTILITY_ALGORITHM_HPP
+#define GAPP_UTILITY_ALGORITHM_HPP
 
+#include "small_vector.hpp"
 #include "type_traits.hpp"
 #include "utility.hpp"
 #include <vector>
@@ -34,9 +35,9 @@ namespace gapp::detail
         return static_cast<CastType>(u);
     }
 
-    inline std::vector<size_t> index_vector(size_t n, size_t first = 0)
+    constexpr small_vector<size_t> index_vector(size_t n, size_t first = 0)
     {
-        std::vector<size_t> indices(n);
+        small_vector<size_t> indices(n);
         std::iota(indices.begin(), indices.end(), first);
 
         return indices;
@@ -44,7 +45,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_reference_t<Iter>, std::iter_reference_t<Iter>>
-    std::vector<size_t> argsort(Iter first, Iter last, Comp&& comp = {})
+    small_vector<size_t> argsort(Iter first, Iter last, Comp&& comp = {})
     {
         GAPP_ASSERT(std::distance(first, last) >= 0);
 
@@ -71,7 +72,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_reference_t<Iter>, std::iter_reference_t<Iter>>
-    std::vector<size_t> partial_argsort(Iter first, Iter middle, Iter last, Comp&& comp = {})
+    small_vector<size_t> partial_argsort(Iter first, Iter middle, Iter last, Comp&& comp = {})
     {
         GAPP_ASSERT(std::distance(first, middle) >= 0);
         GAPP_ASSERT(std::distance(middle, last) >= 0);
@@ -233,7 +234,7 @@ namespace gapp::detail
 
         using ValueType = std::remove_cvref_t<std::iter_reference_t<Iter>>;
 
-        std::vector<ValueType> result;
+        small_vector<ValueType> result;
         result.reserve(last - first);
 
         for (; first != last; ++first)
@@ -265,9 +266,9 @@ namespace gapp::detail
 
     template<std::ranges::random_access_range R, typename Pred>
     requires std::predicate<Pred&, detail::const_reference_t<R>>
-    std::vector<size_t> find_indices(const R& range, Pred&& pred)
+    small_vector<size_t> find_indices(const R& range, Pred&& pred)
     {
-        std::vector<size_t> indices;
+        small_vector<size_t> indices;
 
         for (size_t i = 0; i < range.size(); i++)
         {
@@ -312,29 +313,27 @@ namespace gapp::detail
     }
 
     template<std::ranges::random_access_range R>
-    auto select(R&& range, const std::vector<size_t>& indices)
+    auto select(R&& container, std::span<const size_t> indices)
     {
-        using ValueType = detail::value_t<std::remove_cvref_t<R>>;
-
-        std::vector<ValueType> selected;
+        std::remove_cvref_t<R> selected;
         selected.reserve(indices.size());
 
         for (size_t idx : indices)
         {
-            selected.push_back(detail::forward_like<R>(range[idx]));
+            selected.push_back(detail::forward_like<R>(container[idx]));
         }
 
         return selected;
     }
 
-    template<std::totally_ordered T>
-    constexpr void erase_duplicates(std::vector<T>& container)
+    template<std::ranges::random_access_range R>
+    constexpr void erase_duplicates(R& range)
     {
-        std::sort(container.begin(), container.end());
-        const auto last = std::unique(container.begin(), container.end());
-        container.erase(last, container.end());
+        std::sort(range.begin(), range.end());
+        const auto last = std::unique(range.begin(), range.end());
+        range.erase(last, range.end());
     }
 
 } // namespace gapp::detail
 
-#endif // !GA_UTILITY_ALGORITHM_HPP
\ No newline at end of file
+#endif // !GAPP_UTILITY_ALGORITHM_HPP
\ No newline at end of file
diff --git a/src/utility/cone_tree.cpp b/src/utility/cone_tree.cpp
index 9ea42473..dccb94c7 100644
--- a/src/utility/cone_tree.cpp
+++ b/src/utility/cone_tree.cpp
@@ -15,11 +15,13 @@
 
 namespace gapp::detail
 {
-    using iterator = ConeTree::iterator;
+    using iterator       = ConeTree::iterator;
     using const_iterator = ConeTree::const_iterator;
 
-    using Point = ConeTree::Point;
-    using Node = ConeTree::Node;
+    using Point    = ConeTree::Point;
+    using PointRef = ConeTree::PointRef;
+
+    using Node       = ConeTree::Node;
     using FindResult = ConeTree::FindResult;
 
 
@@ -32,8 +34,7 @@ namespace gapp::detail
         for (const Point& point : points) points_.append_row(point);
 
         nodes_.reserve(4 * points_.size() / MAX_LEAF_ELEMENTS);
-        Node root{ .first = 0, .last = points_.size() };
-        nodes_.push_back(root);
+        nodes_.push_back({ .first = 0, .last = points_.size() });
 
         buildTree();
     }
@@ -64,20 +65,19 @@ namespace gapp::detail
 
         const ptrdiff_t range_len = std::distance(first, last);
 
-        Point center(first->begin(), first->end());
+        Point center(*first);
 
         for (++first; first != last; ++first)
         {
             std::transform(center.begin(), center.end(), first->begin(), center.begin(), std::plus{});
         }
-
         std::transform(center.begin(), center.end(), center.begin(), detail::divide_by(range_len));
 
         return center;
     }
 
     /* 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)
+    static inline double findRadius(const_iterator first, const_iterator last, PointRef center)
     {
         auto distance = std::bind_front(math::euclideanDistanceSq, center);
         auto furthest = detail::max_element(first, last, distance);
@@ -92,7 +92,7 @@ namespace gapp::detail
     }
 
     /* Return the max possible inner product between the point and a point inside the node. */
-    static inline double innerProductUpperBound(const Node& node, const Point& point, double point_norm)
+    static inline double innerProductUpperBound(const Node& node, PointRef point, double point_norm)
     {
         const double center_prod = std::inner_product(point.begin(), point.end(), node.center.begin(), 0.0);
 
@@ -100,7 +100,7 @@ namespace gapp::detail
     }
 
     /* Find the best match in the range [first, last) using linear search. */
-    static FindResult findBestMatchLinear(const Point& query_point, const_iterator first, const_iterator last)
+    static FindResult findBestMatchLinear(PointRef query_point, const_iterator first, const_iterator last)
     {
         GAPP_ASSERT(std::distance(first, last) > 0);
         GAPP_ASSERT(query_point.size() == first->size());
@@ -168,7 +168,7 @@ namespace gapp::detail
         }
     }
 
-    FindResult ConeTree::findBestMatch(const Point& query_point) const
+    FindResult ConeTree::findBestMatch(PointRef query_point) const
     {
         if (points_.empty()) return { points_.end(), 0.0 };
 
diff --git a/src/utility/cone_tree.hpp b/src/utility/cone_tree.hpp
index 27ed7ddf..b2373ed1 100644
--- a/src/utility/cone_tree.hpp
+++ b/src/utility/cone_tree.hpp
@@ -4,6 +4,7 @@
 #define GA_UTILITY_CONE_TREE_HPP
 
 #include "matrix.hpp"
+#include "small_vector.hpp"
 #include <vector>
 #include <span>
 #include <type_traits>
@@ -24,7 +25,8 @@ namespace gapp::detail
         using iterator       = Matrix<double>::iterator;
         using const_iterator = Matrix<double>::const_iterator;
 
-        using Point = std::vector<double>;
+        using Point    = small_vector<double>;
+        using PointRef = std::span<const double>;
 
         struct FindResult
         {
@@ -48,7 +50,7 @@ namespace gapp::detail
         explicit ConeTree(std::span<const Point> points);
 
         /* Returns the closest point in the tree to the query point, and its distance. */
-        FindResult findBestMatch(const Point& query_point) const;
+        FindResult findBestMatch(PointRef query_point) const;
 
         constexpr const_iterator begin() const noexcept { return points_.begin(); }
         constexpr const_iterator end() const noexcept   { return points_.end(); }
diff --git a/src/utility/matrix.hpp b/src/utility/matrix.hpp
index 73636a4a..5a533e2c 100644
--- a/src/utility/matrix.hpp
+++ b/src/utility/matrix.hpp
@@ -3,6 +3,7 @@
 #ifndef GA_UTILITY_MATRIX_HPP
 #define GA_UTILITY_MATRIX_HPP
 
+#include "small_vector.hpp"
 #include "iterators.hpp"
 #include "functional.hpp"
 #include "type_traits.hpp"
@@ -262,6 +263,8 @@ namespace gapp::detail
         template<typename alloc_type>
         operator std::vector<value_type, alloc_type>() const { return std::vector<value_type, alloc_type>(begin(), end()); }
 
+        operator small_vector<value_type>() const { return small_vector<value_type>(begin(), end()); }
+
         operator std::span<copy_const_t<MatrixType, value_type>>() const { return { begin(), end() }; }
 
         constexpr friend bool operator==(const Derived& lhs, const Derived& rhs)
diff --git a/src/utility/qrng.hpp b/src/utility/qrng.hpp
index 1868a081..97ad7387 100644
--- a/src/utility/qrng.hpp
+++ b/src/utility/qrng.hpp
@@ -4,7 +4,7 @@
 #define GA_UTILITY_QRNG_HPP
 
 #include "bounded_value.hpp"
-#include <vector>
+#include "small_vector.hpp"
 #include <concepts>
 #include <cstddef>
 
@@ -18,8 +18,8 @@ namespace gapp::rng
     class QuasiRandom
     {
     public:
-        using result_type = std::vector<RealType>;
-        using state_type  = std::vector<RealType>;
+        using result_type = small_vector<RealType>;
+        using state_type  = small_vector<RealType>;
         using size_type   = std::size_t;
 
         /** Create a quasi-random number generator in @p dim dimensions. */
diff --git a/src/utility/small_vector.hpp b/src/utility/small_vector.hpp
index f2995390..1f551d70 100644
--- a/src/utility/small_vector.hpp
+++ b/src/utility/small_vector.hpp
@@ -5,6 +5,7 @@
 
 #include "scope_exit.hpp"
 #include "utility.hpp"
+#include <vector>
 #include <array>
 #include <algorithm>
 #include <iterator>
@@ -18,7 +19,7 @@
 #include <cstdlib>
 #include <cstddef>
 
-// NOLINTBEGIN(*pointer-arithmetic, *union-access)
+// NOLINTBEGIN(*pointer-arithmetic, *union-access, *redundant-expression, *exception-escape)
 
 namespace gapp::detail
 {
@@ -877,6 +878,8 @@ namespace gapp
         //               OTHER               //
         //-----------------------------------//
 
+        constexpr std::vector<T> std_vec() const { return std::vector<T>{ begin(), end() }; }
+
         constexpr allocator_type get_allocator() const noexcept(std::is_nothrow_copy_constructible_v<A>) { return alloc_; }
 
         constexpr friend bool operator==(const small_vector& lhs, const small_vector& rhs) noexcept
@@ -1076,6 +1079,6 @@ namespace gapp
 
 } // namespace gapp
 
-// NOLINTEND(*pointer-arithmetic, *union-access)
+// NOLINTEND(*pointer-arithmetic, *union-access, *redundant-expression, *exception-escape)
 
 #endif // !GAPP_UTILITY_SMALL_VECTOR_HPP
diff --git a/test/unit/algorithm.cpp b/test/unit/algorithm.cpp
index 3a0fbfdf..c8bb77d2 100644
--- a/test/unit/algorithm.cpp
+++ b/test/unit/algorithm.cpp
@@ -5,6 +5,7 @@
 #include <vector>
 #include <tuple>
 #include <cstddef>
+#include "utility/small_vector.hpp"
 #include "utility/algorithm.hpp"
 #include "utility/rng.hpp"
 #include "utility/utility.hpp"
@@ -63,30 +64,30 @@ TEST_CASE("decrement_mod", "[algorithm]")
 
 TEST_CASE("index_vector", "[algorithm]")
 {
-    REQUIRE(detail::index_vector(3) == std::vector<size_t>{ 0, 1, 2 });
-    REQUIRE(detail::index_vector(4, 2) == std::vector<size_t>{ 2, 3, 4, 5 });
+    REQUIRE(detail::index_vector(3) == small_vector<size_t>{ 0, 1, 2 });
+    REQUIRE(detail::index_vector(4, 2) == small_vector<size_t>{ 2, 3, 4, 5 });
 }
 
 TEST_CASE("argsort", "[algorithm]")
 {
-    const std::vector nums = { 4.0, 0.0, 2.0, 1.0 };
+    const small_vector nums = { 4.0, 0.0, 2.0, 1.0 };
 
     SECTION("iterators")
     {
         const auto indices = detail::argsort(nums.begin(), nums.end());
-        REQUIRE(indices == std::vector<size_t>{ 1, 3, 2, 0 });
+        REQUIRE(indices == small_vector<size_t>{ 1, 3, 2, 0 });
     }
 
     SECTION("reverse iterators")
     {
         const auto indices = detail::argsort(nums.rbegin(), nums.rend());
-        REQUIRE(indices == std::vector<size_t>{ 0, 2, 3, 1 });
+        REQUIRE(indices == small_vector<size_t>{ 0, 2, 3, 1 });
     }
 
     SECTION("compare function")
     {
         const auto indices = detail::argsort(nums.begin(), nums.end(), std::greater<>{});
-        REQUIRE(indices == std::vector<size_t>{ 0, 2, 3, 1 });
+        REQUIRE(indices == small_vector<size_t>{ 0, 2, 3, 1 });
     }
 
     SECTION("empty range")
@@ -98,7 +99,7 @@ TEST_CASE("argsort", "[algorithm]")
 
 TEST_CASE("partial_argsort", "[algorithm]")
 {
-    const std::vector nums = { 4.0, 0.0, 2.0, 1.0, 5.0 };
+    const small_vector nums = { 4.0, 0.0, 2.0, 1.0, 5.0 };
 
     SECTION("iterators")
     {
@@ -136,7 +137,7 @@ TEST_CASE("partial_argsort", "[algorithm]")
 
 TEST_CASE("max_element", "[algorithm]")
 {
-    const std::vector nums = { 4.0, 0.0, 2.0, 5.0, 1.0 };
+    const small_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);
@@ -148,7 +149,7 @@ TEST_CASE("max_element", "[algorithm]")
 
 TEST_CASE("min_element", "[algorithm]")
 {
-    const std::vector nums = { 4.0, 0.0, 2.0, 5.0, 1.0 };
+    const small_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);
@@ -160,7 +161,7 @@ TEST_CASE("min_element", "[algorithm]")
 
 TEST_CASE("argmax", "[algorithm]")
 {
-    const std::vector nums = { 4.0, 0.0, 2.0, 5.0, 1.0 };
+    const small_vector nums = { 4.0, 0.0, 2.0, 5.0, 1.0 };
 
     REQUIRE(detail::argmax(nums.begin(), nums.end()) == 3);
     REQUIRE(detail::argmax(nums.rbegin(), nums.rend()) == 3);
@@ -224,7 +225,7 @@ TEST_CASE("partial_shuffle", "[algorithm]")
 
 TEST_CASE("contains", "[algorithm]")
 {
-    const std::vector nums = { 4.0, 0.0, 2.0, 5.0, 1.0 };
+    const small_vector nums = { 4.0, 0.0, 2.0, 5.0, 1.0 };
 
     REQUIRE(detail::contains(nums.begin(), nums.end(), 0.0));
     REQUIRE(detail::contains(nums.begin(), nums.end(), 1.0));
@@ -235,10 +236,10 @@ TEST_CASE("contains", "[algorithm]")
 
 TEST_CASE("find_all", "[algorithm]")
 {
-    const std::vector nums = { 4, 0, 2, 5, 1 };
+    const small_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{ 5, 1 });
+    REQUIRE(odd_nums == small_vector{ 5, 1 });
 
     const auto big_nums = detail::find_all(nums.begin(), nums.end(), is_big);
     REQUIRE(big_nums.empty());
@@ -249,16 +250,16 @@ TEST_CASE("find_all", "[algorithm]")
 
 TEST_CASE("find_indices", "[algorithm]")
 {
-    const std::vector nums = { 4, 0, 2, 5, 1 };
+    const small_vector nums = { 4, 0, 2, 5, 1 };
 
     const auto odd_num_idxs = detail::find_indices(nums, is_odd);
-    REQUIRE(odd_num_idxs == std::vector{ 3_sz, 4_sz });
+    REQUIRE(odd_num_idxs == small_vector{ 3_sz, 4_sz });
 
     const auto big_num_idxs = detail::find_indices(nums, is_big);
     REQUIRE(big_num_idxs.empty());
 
     const auto all = detail::find_indices(nums, always_true);
-    REQUIRE(all == std::vector{ 0_sz, 1_sz, 2_sz, 3_sz, 4_sz });
+    REQUIRE(all == small_vector{ 0_sz, 1_sz, 2_sz, 3_sz, 4_sz });
 
     const auto none = detail::find_indices(nums, always_false);
     REQUIRE(none.empty());
@@ -266,7 +267,7 @@ TEST_CASE("find_indices", "[algorithm]")
 
 TEST_CASE("index_of", "[algorithm]")
 {
-    const std::vector nums = { 4, 0, 2, 5, 1 };
+    const small_vector nums = { 4, 0, 2, 5, 1 };
 
     REQUIRE(detail::index_of(nums, 4) == 0_sz);
     REQUIRE(detail::index_of(nums, 2) == 2_sz);
@@ -276,7 +277,7 @@ TEST_CASE("index_of", "[algorithm]")
 
 TEST_CASE("find_index", "[algorithm]")
 {
-    const std::vector nums = { 4, 0, 2, 5, 1 };
+    const small_vector nums = { 4, 0, 2, 5, 1 };
 
     const auto first_idx = detail::find_index(nums, always_true);
     REQUIRE(first_idx == 0_sz);
@@ -331,16 +332,16 @@ TEST_CASE("select", "[algorithm]")
 {
     const std::vector nums = { 4, 0, 2, 5, 1, 3, 1 };
 
-    auto selected = detail::select(nums, { 0, 1, 4 });
+    auto selected = detail::select(nums, std::vector{ 0_sz, 1_sz, 4_sz });
     REQUIRE(selected == std::vector{ 4, 0, 1 });
 
-    selected = detail::select(selected, { 2 });
+    selected = detail::select(selected, std::vector{ 2_sz });
     REQUIRE(selected == std::vector{ 1 });
 
     selected = detail::select(nums, { });
     REQUIRE(selected.empty());
 
-    selected = detail::select(std::vector{ 1, 3, 5 }, { 0, 1 });
+    selected = detail::select(std::vector{ 1, 3, 5 }, std::vector{ 0_sz, 1_sz });
     REQUIRE(selected == std::vector{ 1, 3 });
 }
 
diff --git a/test/unit/cone_tree.cpp b/test/unit/cone_tree.cpp
index 2bc2e1b7..1f94049f 100644
--- a/test/unit/cone_tree.cpp
+++ b/test/unit/cone_tree.cpp
@@ -52,10 +52,10 @@ TEST_CASE("cone_tree lookup", "[cone_tree]")
 
     ConeTree tree(points);
 
-    auto best = tree.findBestMatch({ 1.0, 1.0, 0.1 });
+    auto best = tree.findBestMatch(ConeTree::Point{ 1.0, 1.0, 0.1 });
     REQUIRE(*best.elem == ConeTree::Point{ 0.8, 0.8, 0.6 });
 
-    best = tree.findBestMatch({ 0.1, 0.5, 0.8 });
+    best = tree.findBestMatch(ConeTree::Point{ 0.1, 0.5, 0.8 });
     REQUIRE(*best.elem == ConeTree::Point{ 0.6, 0.8, 0.8 });
 }
 
@@ -63,7 +63,7 @@ TEST_CASE("empty_cone_tree", "[cone_tree]")
 {
     ConeTree tree;
 
-    auto best = tree.findBestMatch({ 1.0, 1.0 });
+    auto best = tree.findBestMatch(ConeTree::Point{ 1.0, 1.0 });
 
     REQUIRE(best.elem == tree.end());
 }
\ No newline at end of file
diff --git a/test/unit/pareto_front.cpp b/test/unit/pareto_front.cpp
index 57110ae7..7768425a 100644
--- a/test/unit/pareto_front.cpp
+++ b/test/unit/pareto_front.cpp
@@ -4,6 +4,7 @@
 #include <catch2/catch_template_test_macros.hpp>
 #include <catch2/matchers/catch_matchers_all.hpp>
 #include "core/population.hpp"
+#include "utility/small_vector.hpp"
 #include "utility/utility.hpp"
 #include <vector>
 
@@ -22,7 +23,7 @@ TEST_CASE("find_pareto_front_1D", "[pareto_front]")
 
         auto optimal_indices = findParetoFront1D(fmat);
 
-        REQUIRE(optimal_indices == std::vector{ 1_sz });
+        REQUIRE(optimal_indices == small_vector{ 1_sz });
     }
 
     SECTION("multiple optimum")
@@ -31,7 +32,7 @@ TEST_CASE("find_pareto_front_1D", "[pareto_front]")
 
         auto optimal_indices = findParetoFront1D(fmat);
 
-        REQUIRE(optimal_indices == std::vector{ 2_sz, 7_sz, 10_sz });
+        REQUIRE(optimal_indices == small_vector{ 2_sz, 7_sz, 10_sz });
     }
 
     SECTION("multiple optimum approx")
@@ -40,7 +41,7 @@ TEST_CASE("find_pareto_front_1D", "[pareto_front]")
 
         auto optimal_indices = findParetoFront1D(fmat);
 
-        REQUIRE(optimal_indices == std::vector{ 2_sz, 3_sz, 5_sz, 7_sz, 10_sz });
+        REQUIRE(optimal_indices == small_vector{ 2_sz, 3_sz, 5_sz, 7_sz, 10_sz });
     }
 }
 
@@ -70,14 +71,14 @@ TEMPLATE_TEST_CASE_SIG("find_pareto_front_nd", "[pareto_front]", ((auto F), F),
 
         auto optimal_indices = F(fmat);
 
-        REQUIRE(optimal_indices == std::vector<size_t>{ 3 });
+        REQUIRE(optimal_indices == small_vector<size_t>{ 3 });
     }
 
     SECTION("multiple optimum")
     {
         ScopedTolerances _(0.0, 0.0);
 
-        auto optimal_indices = F(fmat);
+        auto optimal_indices = F(fmat).std_vec();
 
         REQUIRE_THAT(optimal_indices, Matchers::UnorderedEquals(std::vector<size_t>{ 5, 9, 11, 13, 14 }));
     }
@@ -86,7 +87,7 @@ TEMPLATE_TEST_CASE_SIG("find_pareto_front_nd", "[pareto_front]", ((auto F), F),
     {
         ScopedTolerances _(0.1, 0.0);
 
-        auto optimal_indices = F(fmat);
+        auto optimal_indices = F(fmat).std_vec();
 
         REQUIRE_THAT(optimal_indices, Matchers::UnorderedEquals(std::vector<size_t>{ 4, 5, 9, 11, 12, 13, 14 }));
     }
diff --git a/test/unit/problems.cpp b/test/unit/problems.cpp
index 497b7f24..efae98da 100644
--- a/test/unit/problems.cpp
+++ b/test/unit/problems.cpp
@@ -42,7 +42,7 @@ TEMPLATE_TEST_CASE("single_objective_problems", "[problems]", Sphere, Rastrigin,
 
     TestType func(var_count);
 
-    REQUIRE_THAT( func(func.optimum()), Approx(func.optimal_value()).margin(1E-6) );
+    REQUIRE_THAT( func(func.optimum()).std_vec(), Approx(func.optimal_value().std_vec()).margin(1E-6) );
 
     const auto random_chrom = randomChromosome(func.bounds());
 
@@ -58,7 +58,7 @@ TEST_CASE("kursawe", "[problems]")
 
     Kursawe func(var_count);
 
-    REQUIRE_THAT( func(func.optimum()), Approx(func.optimal_value()).margin(1E-6) );
+    REQUIRE_THAT( func(func.optimum()).std_vec(), Approx(func.optimal_value().std_vec()).margin(1E-6));
 
     REQUIRE( !paretoCompareLess(func.ideal_point(), func.nadir_point()) );
     REQUIRE( !paretoCompareLess(func.optimal_value(), func.nadir_point()) );
@@ -80,7 +80,7 @@ TEMPLATE_TEST_CASE("zdt_suite", "[problems]", ZDT1, ZDT2, ZDT3, ZDT4, ZDT5, ZDT6
 
     TestType func(var_count);
 
-    REQUIRE_THAT(func(func.optimum()), Approx(func.optimal_value()).margin(1E-6));
+    REQUIRE_THAT(func(func.optimum()).std_vec(), Approx(func.optimal_value().std_vec()).margin(1E-6));
 
     REQUIRE( !paretoCompareLess(func.ideal_point(), func.nadir_point()) );
     REQUIRE( !paretoCompareLess(func.optimal_value(), func.nadir_point()) );
@@ -102,7 +102,7 @@ TEMPLATE_TEST_CASE("dtlz_suite", "[problems]", DTLZ1, DTLZ2, DTLZ3, DTLZ4, DTLZ5
 
     TestType func(num_obj);
 
-    REQUIRE_THAT( func(func.optimum()), Approx(func.optimal_value()).margin(1E-6) );
+    REQUIRE_THAT( func(func.optimum()).std_vec(), Approx(func.optimal_value().std_vec()).margin(1E-6) );
 
     REQUIRE( !paretoCompareLess(func.ideal_point(), func.nadir_point()) );
     REQUIRE( !paretoCompareLess(func.optimal_value(), func.nadir_point()) );
diff --git a/test/unit/replacement.cpp b/test/unit/replacement.cpp
index 2a0945ae..1ea18d20 100644
--- a/test/unit/replacement.cpp
+++ b/test/unit/replacement.cpp
@@ -44,7 +44,7 @@ TEST_CASE("replacement_best", "[replacement][single-objective]")
     math::ScopedTolerances _(0.0, 0.0);
 
     std::unique_ptr<Replacement> replacement = std::make_unique<KeepBest>();
-    const auto indices = replacement->nextPopulationImpl(context, fitness_matrix);
+    const auto indices = replacement->nextPopulationImpl(context, fitness_matrix).std_vec();
 
     const std::vector<size_t> expected = { 0, 1, 4, 5, 9, 11, 12, 15, 16, 17 };
 
@@ -54,7 +54,7 @@ TEST_CASE("replacement_best", "[replacement][single-objective]")
 TEST_CASE("replacement_children", "[replacement][single-objective]")
 {
     std::unique_ptr<Replacement> replacement = std::make_unique<KeepChildren>();
-    const auto indices = replacement->nextPopulationImpl(context, fitness_matrix);
+    const auto indices = replacement->nextPopulationImpl(context, fitness_matrix).std_vec();
 
     const std::vector<size_t> expected = { 10, 11, 12, 13, 14, 15, 16, 17, 18, 19 };
 
@@ -64,7 +64,7 @@ TEST_CASE("replacement_children", "[replacement][single-objective]")
 TEST_CASE("replacement_elitism", "[replacement][single-objective]")
 {
     std::unique_ptr<Replacement> replacement = std::make_unique<Elitism>(2);
-    const auto indices = replacement->nextPopulationImpl(context, fitness_matrix);
+    const auto indices = replacement->nextPopulationImpl(context, fitness_matrix).std_vec();
 
     const std::vector<size_t> expected = { 0, 5, 10, 11, 12, 13, 14, 15, 16, 17 };