Possible 50% performance improvement for NSGA-II on MSVC

include/pagmo/detail/custom_comparisons.hpp

@@ -54,7 +54,38 @@ namespace detail
 template <typename T, bool After = true>
 inline bool less_than_f(T a, T b)
 {
-    static_assert(std::is_floating_point<T>::value, "less_than_f can be used only with floating-point types.");
+    static_assert(std::is_floating_point_v<T>, "less_than_f can be used only with floating-point types.");
+
+#if defined(__cpp_impl_three_way_comparison)
+
+    const auto cmp = a <=> b;
+    // True or false returned immediately, if there's no NaN (std::isnan is very slow)
+    if (cmp == std::partial_ordering::less) {
+        return true;
+    }
+    if (cmp == std::partial_ordering::greater || cmp == std::partial_ordering::equivalent) {
+        return false;
+    }
+
+    // Result was std::partial_ordering::unordered, one of the operands must be NaN
+    const bool a_nan = std::isnan(a);
+    const bool b_nan = std::isnan(b);
+
+    // If only B is NaN, After must be returned
+    if (!a_nan && b_nan) {
+        return After;
+    }
+    // If only A is NaN, !After must be returned
+    if (a_nan && !b_nan) {
+        return !After;
+    }
+
+    // If both are NaN, false must be returned
+    return false;
+
+#else
+
+    // Pre-C++20 implementation
     if (!std::isnan(a)) {
         if (!std::isnan(b))
             return a < b; // a < b
@@ -66,6 +97,8 @@ inline bool less_than_f(T a, T b)
         else
             return false; // nan < nan
     }
+
+#endif
 }
 
 // Greater than compares floating point types placing nans after inf or before -inf
@@ -74,7 +107,33 @@ inline bool less_than_f(T a, T b)
 template <typename T, bool After = true>
 inline bool greater_than_f(T a, T b)
 {
-    static_assert(std::is_floating_point<T>::value, "greater_than_f can be used only with floating-point types.");
+    static_assert(std::is_floating_point_v<T>, "greater_than_f can be used only with floating-point types.");
+
+#if defined(__cpp_impl_three_way_comparison)
+
+    auto cmp = a <=> b;
+    if (cmp == std::partial_ordering::greater) {
+        return true;
+    }
+    if (cmp == std::partial_ordering::less || cmp == std::partial_ordering::equivalent) {
+        return false;
+    }
+
+    const bool a_nan = std::isnan(a);
+    const bool b_nan = std::isnan(b);
+
+    if (!a_nan && b_nan) {
+        return !After;
+    }
+    if (a_nan && !b_nan) {
+        return After;
+    }
+
+    return false;
+
+#else
+
+    // Pre-C++20 implementation
     if (!std::isnan(a)) {
         if (!std::isnan(b))
             return a > b; // a > b
@@ -86,17 +145,38 @@ inline bool greater_than_f(T a, T b)
         else
             return false; // nan > nan
     }
+
+#endif
 }
 
 // equal_to than compares floating point types considering nan==nan
 template <typename T>
 inline bool equal_to_f(T a, T b)
 {
+#if defined(__cpp_impl_three_way_comparison)
+
+    auto cmp = a <=> b;
+    if (cmp == std::partial_ordering::equivalent) {
+        return true;
+    }
+
+    if (cmp != std::partial_ordering::unordered) {
+        return false;
+    }
+
+    return std::isnan(a) && std::isnan(b);
+
+#else
+
+    // Pre-C++20 implementation
     static_assert(std::is_floating_point<T>::value, "equal_to_f can be used only with floating-point types.");
     if (!std::isnan(a) && !std::isnan(b)) {
         return a == b;
     }
     return std::isnan(a) && std::isnan(b);
+
+#endif
+
 }
 
 // equal_to_vf than compares vectors of floating point types considering nan==nan

include/pagmo/utils/multi_objective.hpp

@@ -59,7 +59,7 @@ PAGMO_DLL_PUBLIC void reksum(std::vector<std::vector<double>> &, const std::vect
 } // namespace detail
 
 // Pareto-dominance
-PAGMO_DLL_PUBLIC bool pareto_dominance(const vector_double &, const vector_double &);
+PAGMO_DLL_PUBLIC inline bool pareto_dominance(const vector_double &, const vector_double &);
 
 // Non dominated front 2D (Kung's algorithm)
 PAGMO_DLL_PUBLIC std::vector<pop_size_t> non_dominated_front_2d(const std::vector<vector_double> &);
@@ -71,6 +71,9 @@ using fnds_return_type = std::tuple<std::vector<std::vector<pop_size_t>>, std::v
 // Fast non dominated sorting
 PAGMO_DLL_PUBLIC fnds_return_type fast_non_dominated_sorting(const std::vector<vector_double> &);
 
+// Fast non dominated sorting with pre-allocated buffers
+PAGMO_DLL_PUBLIC void fast_non_dominated_sorting_buffered(const std::vector<vector_double> &, fnds_return_type&);
+
 // Crowding distance
 PAGMO_DLL_PUBLIC vector_double crowding_distance(const std::vector<vector_double> &);
 
@@ -80,6 +83,9 @@ PAGMO_DLL_PUBLIC std::vector<pop_size_t> sort_population_mo(const std::vector<ve
 // Selects the best N individuals in multi-objective optimization
 PAGMO_DLL_PUBLIC std::vector<pop_size_t> select_best_N_mo(const std::vector<vector_double> &, pop_size_t);
 
+// Selects the best N individuals in multi-objective optimization with pre-allocated buffers
+PAGMO_DLL_PUBLIC std::vector<pop_size_t> select_best_N_mo_buffered(const std::vector<vector_double> &, pop_size_t, fnds_return_type&);
+
 // Ideal point
 PAGMO_DLL_PUBLIC vector_double ideal(const std::vector<vector_double> &);
 

src/algorithms/nsga2.cpp

@@ -136,6 +136,10 @@ population nsga2::evolve(population pop) const
     vector_double::size_type parent1_idx, parent2_idx;
     std::pair<vector_double, vector_double> children;
 
+    // We're using select_best_N_mo_buffered(), which prevents re-allocation of this structure
+    fnds_return_type best_N_buffer{};
+    fnds_return_type fnds_buffer{};
+
     std::iota(shuffle1.begin(), shuffle1.end(), vector_double::size_type(0));
     std::iota(shuffle2.begin(), shuffle2.end(), vector_double::size_type(0));
 
@@ -181,10 +185,10 @@ population nsga2::evolve(population pop) const
         std::shuffle(shuffle2.begin(), shuffle2.end(), m_e);
 
         // 1 - We compute crowding distance and non dominated rank for the current population
-        auto fnds_res = fast_non_dominated_sorting(pop.get_f());
-        auto ndf = std::get<0>(fnds_res); // non dominated fronts [[0,3,2],[1,5,6],[4],...]
+        fast_non_dominated_sorting_buffered(pop.get_f(), fnds_buffer);
+        auto& ndf = std::get<0>(fnds_buffer); // non dominated fronts [[0,3,2],[1,5,6],[4],...]
         vector_double pop_cd(NP);         // crowding distances of the whole population
-        auto ndr = std::get<3>(fnds_res); // non domination rank [0,1,0,0,2,1,1, ... ]
+        auto& ndr = std::get<3>(fnds_buffer); // non domination rank [0,1,0,0,2,1,1, ... ]
         for (const auto &front_idxs : ndf) {
             if (front_idxs.size() == 1u) { // handles the case where the front has collapsed to one point
                 pop_cd[front_idxs[0]] = std::numeric_limits<double>::infinity();
@@ -297,7 +301,7 @@ population nsga2::evolve(population pop) const
         }
         // This method returns the sorted N best individuals in the population according to the crowded comparison
         // operator
-        best_idx = select_best_N_mo(popnew.get_f(), NP);
+        best_idx = select_best_N_mo_buffered(popnew.get_f(), NP, best_N_buffer);
         // We insert into the population
         for (population::size_type i = 0; i < NP; ++i) {
             pop.set_xf(i, popnew.get_x()[best_idx[i]], popnew.get_f()[best_idx[i]]);

src/utils/multi_objective.cpp

@@ -94,7 +94,7 @@ void reksum(std::vector<std::vector<double>> &retval, const std::vector<pop_size
  *
  * @throws std::invalid_argument if the dimensions of the two objectives are different
  */
-bool pareto_dominance(const vector_double &obj1, const vector_double &obj2)
+inline bool pareto_dominance(const vector_double &obj1, const vector_double &obj2)
 {
     if (obj1.size() != obj2.size()) {
         pagmo_throw(std::invalid_argument,
@@ -256,6 +256,84 @@ fnds_return_type fast_non_dominated_sorting(const std::vector<vector_double> &po
                            std::move(non_dom_rank));
 }
 
+/// Fast non dominated sorting with pre-allocated buffers
+/**
+ * Improved version of fast_non_dominated_sorting(), doesn't allocate any vectors inside the function body, instead lets
+ * user pass an object of type fnds_return_type, as one of the parameters. This improves performance when called inside
+ * hot path loops.
+ *
+ * @param points An std::vector containing the objectives of different individuals. Example
+ * {{1,2,3},{-2,3,7},{-1,-2,-3},{0,0,0}}
+ * @param sorting_results Return values of this function, same as fast_non_dominated_sorting()
+ *
+ * @throws std::invalid_argument If the size of \p points is not at least 2
+ */
+void fast_non_dominated_sorting_buffered(const std::vector<vector_double> &points, fnds_return_type& sorting_results)
+{
+    auto N = points.size();
+    // We make sure to have two points at least (one could also be allowed)
+    if (N < 2u) {
+        pagmo_throw(std::invalid_argument, "At least two points are needed for fast_non_dominated_sorting: "
+                                               + std::to_string(N) + " detected.");
+    }
+
+    // These are the return values, we need to invalidate previous content of these vectors
+    std::vector<std::vector<pop_size_t>>& non_dom_fronts = std::get<0>(sorting_results);
+    std::vector<std::vector<pop_size_t>>& dom_list = std::get<1>(sorting_results);
+    std::vector<pop_size_t>& dom_count = std::get<2>(sorting_results);
+    std::vector<pop_size_t>& non_dom_rank = std::get<3>(sorting_results);
+    non_dom_fronts.clear();
+    non_dom_fronts.resize(1u);
+
+    dom_count.clear();
+    dom_list.resize(N);
+
+    dom_count.assign(N, 0u);
+    non_dom_rank.assign(N, 0u);
+
+    // Start the fast non dominated sort algorithm
+    for (decltype(N) i = 0u; i < N; ++i) {
+        dom_list[i].clear();
+        dom_count[i] = 0u;
+        for (decltype(N) j = 0u; j < i; ++j) {
+            if (pareto_dominance(points[i], points[j])) {
+                dom_list[i].push_back(j);
+                ++dom_count[j];
+            } else if (pareto_dominance(points[j], points[i])) {
+                dom_list[j].push_back(i);
+                ++dom_count[i];
+            }
+        }
+    }
+    for (decltype(N) i = 0u; i < N; ++i) {
+        if (dom_count[i] == 0u) {
+            non_dom_rank[i] = 0u;
+            non_dom_fronts[0].push_back(i);
+        }
+    }
+    // we copy dom_count as we want to output its value at this point
+    auto dom_count_copy(dom_count);
+    auto current_front = non_dom_fronts[0];
+    std::vector<std::vector<pop_size_t>>::size_type front_counter(0u);
+    while (!current_front.empty()) {
+        std::vector<pop_size_t> next_front;
+        for (const decltype(current_front.size()) p : current_front) {
+            for (const decltype(dom_list[p].size()) q : dom_list[p]) {
+                --dom_count_copy[q];
+                if (dom_count_copy[q] == 0u) {
+                    non_dom_rank[q] = front_counter + 1u;
+                    next_front.push_back(q);
+                }
+            }
+        }
+        ++front_counter;
+        current_front = next_front;
+        if (!current_front.empty()) {
+            non_dom_fronts.push_back(current_front);
+        }
+    }
+}
+
 /// Crowding distance
 /**
  * An implementation of the crowding distance. Complexity is \f$ O(MNlog(N))\f$ where \f$M\f$ is the number of
@@ -395,6 +473,65 @@ std::vector<pop_size_t> select_best_N_mo(const std::vector<vector_double> &input
     return retval;
 }
 
+/// Selects the best N individuals in multi-objective optimization
+/**
+ * select_best_N_mo() with pre-allocated buffers.
+ */
+std::vector<pop_size_t> select_best_N_mo_buffered(const std::vector<vector_double> &input_f, pop_size_t N, fnds_return_type& fnds_buffer)
+{
+    if (N == 0u) { // corner case
+        return {};
+    }
+    if (input_f.size() == 0u) { // corner case
+        return {};
+    }
+    if (input_f.size() == 1u) { // corner case
+        return {0u};
+    }
+    if (N >= input_f.size()) { // corner case
+        std::vector<pop_size_t> retval(input_f.size());
+        std::iota(retval.begin(), retval.end(), pop_size_t(0u));
+        return retval;
+    }
+    std::vector<pop_size_t> retval;
+    std::vector<pop_size_t>::size_type front_id(0u);
+    // Run fast-non-dominated sorting
+    fast_non_dominated_sorting_buffered(input_f, fnds_buffer);
+    // Insert all non dominated fronts if not more than N
+    for (const auto &front : std::get<0>(fnds_buffer)) {
+        if (retval.size() + front.size() <= N) {
+            for (auto i : front) {
+                retval.push_back(i);
+            }
+            if (retval.size() == N) {
+                return retval;
+            }
+            ++front_id;
+        } else {
+            break;
+        }
+    }
+    const auto& front = std::get<0>(fnds_buffer)[front_id];
+    std::vector<vector_double> non_dom_fits(front.size());
+    // Run crowding distance for the front
+    for (decltype(front.size()) i = 0u; i < front.size(); ++i) {
+        non_dom_fits[i] = input_f[front[i]];
+    }
+    vector_double cds(crowding_distance(non_dom_fits));
+    // We now have front and crowding distance, we sort the front w.r.t. the crowding
+    std::vector<pop_size_t> idxs(front.size());
+    std::iota(idxs.begin(), idxs.end(), pop_size_t(0u));
+    std::sort(idxs.begin(), idxs.end(), [&cds](pop_size_t idx1, pop_size_t idx2) {
+        return detail::greater_than_f(cds[idx1], cds[idx2]);
+    }); // Descending order1
+    auto remaining = N - retval.size();
+    for (decltype(remaining) i = 0u; i < remaining; ++i) {
+        retval.push_back(front[idxs[i]]);
+    }
+    return retval;
+}
+
+
 /// Sorts a population in multi-objective optimization
 /**
  * Sorts a population (intended here as an <tt>std::vector<vector_double></tt> containing the  objective vectors)