diff --git a/include/pagmo/detail/custom_comparisons.hpp b/include/pagmo/detail/custom_comparisons.hpp index 58aff766c..156395432 100644 --- a/include/pagmo/detail/custom_comparisons.hpp +++ b/include/pagmo/detail/custom_comparisons.hpp @@ -54,7 +54,38 @@ namespace detail template inline bool less_than_f(T a, T b) { - static_assert(std::is_floating_point::value, "less_than_f can be used only with floating-point types."); + static_assert(std::is_floating_point_v, "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 inline bool greater_than_f(T a, T b) { - static_assert(std::is_floating_point::value, "greater_than_f can be used only with floating-point types."); + static_assert(std::is_floating_point_v, "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 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::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 diff --git a/include/pagmo/utils/multi_objective.hpp b/include/pagmo/utils/multi_objective.hpp index 4c298f727..515dcac1b 100644 --- a/include/pagmo/utils/multi_objective.hpp +++ b/include/pagmo/utils/multi_objective.hpp @@ -59,7 +59,7 @@ PAGMO_DLL_PUBLIC void reksum(std::vector> &, 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 non_dominated_front_2d(const std::vector &); @@ -71,6 +71,9 @@ using fnds_return_type = std::tuple>, std::v // Fast non dominated sorting PAGMO_DLL_PUBLIC fnds_return_type fast_non_dominated_sorting(const std::vector &); +// Fast non dominated sorting with pre-allocated buffers +PAGMO_DLL_PUBLIC void fast_non_dominated_sorting_buffered(const std::vector &, fnds_return_type&); + // Crowding distance PAGMO_DLL_PUBLIC vector_double crowding_distance(const std::vector &); @@ -80,6 +83,9 @@ PAGMO_DLL_PUBLIC std::vector sort_population_mo(const std::vector select_best_N_mo(const std::vector &, pop_size_t); +// Selects the best N individuals in multi-objective optimization with pre-allocated buffers +PAGMO_DLL_PUBLIC std::vector select_best_N_mo_buffered(const std::vector &, pop_size_t, fnds_return_type&); + // Ideal point PAGMO_DLL_PUBLIC vector_double ideal(const std::vector &); diff --git a/src/algorithms/nsga2.cpp b/src/algorithms/nsga2.cpp index 4ddd6b198..2ec0d0f2a 100644 --- a/src/algorithms/nsga2.cpp +++ b/src/algorithms/nsga2.cpp @@ -136,6 +136,10 @@ population nsga2::evolve(population pop) const vector_double::size_type parent1_idx, parent2_idx; std::pair 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::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]]); diff --git a/src/utils/multi_objective.cpp b/src/utils/multi_objective.cpp index 585df9b42..2594214ba 100644 --- a/src/utils/multi_objective.cpp +++ b/src/utils/multi_objective.cpp @@ -94,7 +94,7 @@ void reksum(std::vector> &retval, const std::vector &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 &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>& non_dom_fronts = std::get<0>(sorting_results); + std::vector>& dom_list = std::get<1>(sorting_results); + std::vector& dom_count = std::get<2>(sorting_results); + std::vector& 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>::size_type front_counter(0u); + while (!current_front.empty()) { + std::vector 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 select_best_N_mo(const std::vector &input return retval; } +/// Selects the best N individuals in multi-objective optimization +/** + * select_best_N_mo() with pre-allocated buffers. + */ +std::vector select_best_N_mo_buffered(const std::vector &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 retval(input_f.size()); + std::iota(retval.begin(), retval.end(), pop_size_t(0u)); + return retval; + } + std::vector retval; + std::vector::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 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 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 std::vector containing the objective vectors)