diff --git a/include/pagmo/algorithms/nsga3.hpp b/include/pagmo/algorithms/nsga3.hpp index 225aa2ec8..a6787ec21 100644 --- a/include/pagmo/algorithms/nsga3.hpp +++ b/include/pagmo/algorithms/nsga3.hpp @@ -38,23 +38,23 @@ class PAGMO_DLL_PUBLIC nsga3{ const log_type &get_log() const { return m_log; } void set_verbosity(unsigned level) { m_verbosity = level; } unsigned get_verbosity() const { return m_verbosity; } - void set_seed(unsigned seed) { reng.seed(seed); seed = seed; } - unsigned get_seed() const { return seed; } - bool has_memory() const {return use_memory; } + void set_seed(unsigned seed) { m_reng.seed(seed); m_seed = seed; } + unsigned get_seed() const { return m_seed; } + bool has_memory() const {return m_use_memory; } private: - unsigned ngen; - double cr; // crossover - double eta_c; // eta crossover - double mut; // mutation - double eta_mut; // eta mutation - size_t divisions; // Reference Point hyperplane subdivisions - unsigned seed; // Seed for PRNG initialisation - bool use_memory; // Preserve extremes, ideal, nadir across generations - mutable NSGA3Memory memory{}; - mutable detail::random_engine_type reng; // Defaults to std::mt19937 + unsigned m_gen; + double m_cr; // crossover + double m_eta_c; // eta crossover + double m_mut; // mutation + double m_eta_mut; // eta mutation + size_t m_divisions; // Reference Point hyperplane subdivisions + unsigned m_seed; // Seed for PRNG initialisation + bool m_use_memory; // Preserve extremes, ideal, nadir across generations + mutable NSGA3Memory m_memory{}; + mutable detail::random_engine_type m_reng; // Defaults to std::mt19937 mutable log_type m_log; unsigned m_verbosity; - mutable std::vector refpoints; + mutable std::vector m_refpoints; // Serialisation support friend class boost::serialization::access; template diff --git a/include/pagmo/utils/multi_objective.hpp b/include/pagmo/utils/multi_objective.hpp index 5b8205162..4bd32ea9f 100644 --- a/include/pagmo/utils/multi_objective.hpp +++ b/include/pagmo/utils/multi_objective.hpp @@ -227,7 +227,6 @@ PAGMO_DLL_PUBLIC double perpendicular_distance(const std::vector &, cons /* Choose single random element from vector container */ template -PAGMO_DLL_PUBLIC T choose_random_element(const std::vector &container){ std::vector choice; std::sample(container.begin(), container.end(), std::back_inserter(choice), 1, std::mt19937{random_device::next()}); diff --git a/include/pagmo/utils/reference_point.hpp b/include/pagmo/utils/reference_point.hpp index 422cfb8b4..3cc77fb56 100644 --- a/include/pagmo/utils/reference_point.hpp +++ b/include/pagmo/utils/reference_point.hpp @@ -21,9 +21,9 @@ namespace pagmo{ class PAGMO_DLL_PUBLIC ReferencePoint{ public: ReferencePoint(size_t nobj); - ~ReferencePoint(); size_t dim() const; - double& operator[](int); + //double& operator[](int); + double& operator[](size_t); friend PAGMO_DLL_PUBLIC std::ostream& operator<<(std::ostream& ostr, const ReferencePoint& rp); void increment_members(){ ++nmembers; } void decrement_members(){ --nmembers; } @@ -56,7 +56,7 @@ void associate_with_reference_points( size_t identify_niche_point(std::vector &); -size_t n_choose_k(unsigned, unsigned); +size_t n_choose_k(size_t, size_t); } // namespace pagmo diff --git a/src/algorithms/nsga3.cpp b/src/algorithms/nsga3.cpp index c18d6d3ec..5bd344744 100644 --- a/src/algorithms/nsga3.cpp +++ b/src/algorithms/nsga3.cpp @@ -28,8 +28,8 @@ namespace pagmo{ nsga3::nsga3(unsigned gen, double cr, double eta_c, double mut, double eta_mut, size_t divisions, unsigned seed, bool use_memory) - : ngen(gen), cr(cr), eta_c(eta_c), mut(mut), eta_mut(eta_mut), - divisions(divisions), seed(seed), use_memory(use_memory), reng(seed){ + : m_gen(gen), m_cr(cr), m_eta_c(eta_c), m_mut(mut), m_eta_mut(eta_mut), + m_divisions(divisions), m_seed(seed), m_use_memory(use_memory), m_reng(seed){ // Validate ctor args if(cr < 0.0 || cr > 1.0){ pagmo_throw(std::invalid_argument, "The crossover probability must be in the range [0, 1], while a value of " @@ -63,11 +63,11 @@ nsga3::nsga3(unsigned gen, double cr, double eta_c, double mut, double eta_mut, std::vector nsga3::generate_uniform_reference_points(size_t nobjs, size_t divisions) const{ ReferencePoint rp(nobjs); - if(!refpoints.empty()){ - refpoints.clear(); + if(!m_refpoints.empty()){ + m_refpoints.clear(); } - refpoints = generate_reference_point_level(rp, divisions, 0, divisions); - return refpoints; + m_refpoints = generate_reference_point_level(rp, divisions, 0, divisions); + return m_refpoints; } @@ -77,14 +77,14 @@ std::vector> nsga3::translate_objectives(population pop) con auto objs = pop.get_f(); std::vector p_ideal; if(has_memory()){ - decltype(objs) combined{memory.v_ideal}; - if(memory.v_ideal.size() != 0){ // i.e. not first gen + decltype(objs) combined{m_memory.v_ideal}; + if(m_memory.v_ideal.size() != 0){ // i.e. not first gen combined.insert(combined.end(), objs.begin(), objs.end()); p_ideal = ideal(combined); }else{ p_ideal = ideal(objs); } - memory.v_ideal = p_ideal; + m_memory.v_ideal = p_ideal; }else{ p_ideal = ideal(objs); } @@ -113,16 +113,16 @@ std::vector> nsga3::find_extreme_points(population pop, std::vector min_obj{}; if(has_memory()){ - if(memory.v_extreme.size() == 0){ + if(m_memory.v_extreme.size() == 0){ for(size_t idx=0; idx(nobj, {})); + m_memory.v_extreme.push_back(std::vector(nobj, {})); } }else{ - for(size_t p=0; p> nsga3::find_extreme_points(population pop, } points.push_back(std::vector(min_obj)); if(has_memory()){ - memory.v_extreme[i] = std::vector(min_obj); + m_memory.v_extreme[i] = std::vector(min_obj); } } @@ -195,14 +195,14 @@ std::vector nsga3::find_intercepts(population pop, std::vector v_nadir; if(has_memory()){ - decltype(objs) combined{memory.v_nadir}; - if(memory.v_nadir.size() != 0){ + decltype(objs) combined{m_memory.v_nadir}; + if(m_memory.v_nadir.size() != 0){ combined.insert(combined.end(), objs.begin(), objs.end()); v_nadir = nadir(combined); }else{ v_nadir = nadir(objs); } - memory.v_nadir = v_nadir; + m_memory.v_nadir = v_nadir; }else{ v_nadir = nadir(objs); } @@ -241,8 +241,6 @@ population nsga3::evolve(population pop) const{ auto dim_i = prob.get_nix(); auto NP = pop.size(); - // Initialize the population - /* Verify problem characteristics: * - Has multiple objectives * - Is not stochastic @@ -271,11 +269,11 @@ population nsga3::evolve(population pop) const{ "NSGA-III requires a population greater than 5 and which is divisible by 4." "Detected input population size is: " + std::to_string(NP)); } - size_t num_rps = n_choose_k(prob.get_nf() + divisions - 1, divisions); + size_t num_rps = n_choose_k(prob.get_nf() + m_divisions - 1, m_divisions); if(NP <= num_rps){ pagmo_throw(std::invalid_argument, "Population size must exceed number of reference points. NP = " - + std::to_string(NP) + " while " + std::to_string(divisions) + " divisions for " + + std::to_string(NP) + " while " + std::to_string(m_divisions) + " divisions for " "reference points gives a total of " + std::to_string(num_rps) + " points."); } @@ -289,13 +287,13 @@ population nsga3::evolve(population pop) const{ std::iota(shuffle1.begin(), shuffle1.end(), vector_double::size_type(0)); std::iota(shuffle2.begin(), shuffle2.end(), vector_double::size_type(0)); - for(decltype(ngen)gen = 1u; gen <= ngen; gen++){ + for(decltype(m_gen)gen = 1u; gen <= m_gen; gen++){ // Copy existing population population popnew(pop); // Permute population indices - std::shuffle(shuffle1.begin(), shuffle1.end(), reng); - std::shuffle(shuffle2.begin(), shuffle2.end(), reng); + std::shuffle(shuffle1.begin(), shuffle1.end(), m_reng); + std::shuffle(shuffle2.begin(), shuffle2.end(), m_reng); /* 1. Generate offspring population Q_t * 2. R = P_t U Q_t @@ -304,7 +302,7 @@ population nsga3::evolve(population pop) const{ if(m_verbosity > 0u){ std::vector p_ideal = ideal(pop.get_f()); - if (gen % m_verbosity == 1u || m_verbosity == 1u) { + if (m_gen % m_verbosity == 1u || m_verbosity == 1u) { // We compute the ideal point // Every 50 lines print the column names if (count % 50u == 1u) { @@ -318,7 +316,7 @@ population nsga3::evolve(population pop) const{ } print('\n'); } - print(std::setw(7), gen, std::setw(15), prob.get_fevals() - fevals0); + print(std::setw(7), m_gen, std::setw(15), prob.get_fevals() - fevals0); for (decltype(p_ideal.size()) i = 0u; i < p_ideal.size(); ++i) { if (i >= 5u) { break; @@ -328,18 +326,18 @@ population nsga3::evolve(population pop) const{ print('\n'); ++count; } - m_log.emplace_back(gen, prob.get_fevals() - fevals0, p_ideal); + m_log.emplace_back(m_gen, prob.get_fevals() - fevals0, p_ideal); } // Offspring generation for (decltype(NP) i = 0; i < NP; i += 4) { // We create two offsprings using the shuffled list 1 decltype(shuffle1) parents1; - std::sample(shuffle1.begin(), shuffle1.end(), std::back_inserter(parents1), 2, std::mt19937{reng()}); + std::sample(shuffle1.begin(), shuffle1.end(), std::back_inserter(parents1), 2, std::mt19937{m_reng()}); children = detail::sbx_crossover_impl(pop.get_x()[parents1[0]], pop.get_x()[parents1[1]], bounds, dim_i, - cr, eta_c, reng); - detail::polynomial_mutation_impl(children.first, bounds, dim_i, mut, eta_mut, reng); - detail::polynomial_mutation_impl(children.second, bounds, dim_i, mut, eta_mut, reng); + m_cr, m_eta_c, m_reng); + detail::polynomial_mutation_impl(children.first, bounds, dim_i, m_mut, m_eta_mut, m_reng); + detail::polynomial_mutation_impl(children.second, bounds, dim_i, m_mut, m_eta_mut, m_reng); // Evaluation via prob ensures feval counter is correctly updated auto f1 = prob.fitness(children.first); auto f2 = prob.fitness(children.second); @@ -348,11 +346,11 @@ population nsga3::evolve(population pop) const{ // Repeat with the shuffled list 2 decltype(shuffle2) parents2; - std::sample(shuffle2.begin(), shuffle2.end(), std::back_inserter(parents2), 2, std::mt19937{reng()}); + std::sample(shuffle2.begin(), shuffle2.end(), std::back_inserter(parents2), 2, std::mt19937{m_reng()}); children = detail::sbx_crossover_impl(pop.get_x()[parents2[0]], pop.get_x()[parents2[1]], bounds, dim_i, - cr, eta_c, reng); - detail::polynomial_mutation_impl(children.first, bounds, dim_i, mut, eta_mut, reng); - detail::polynomial_mutation_impl(children.second, bounds, dim_i, mut, eta_mut, reng); + m_cr, m_eta_c, m_reng); + detail::polynomial_mutation_impl(children.first, bounds, dim_i, m_mut, m_eta_mut, m_reng); + detail::polynomial_mutation_impl(children.second, bounds, dim_i, m_mut, m_eta_mut, m_reng); f1 = prob.fitness(children.first); f2 = prob.fitness(children.second); popnew.push_back(children.first, f1); @@ -389,7 +387,7 @@ std::vector nsga3::selection(population &R, size_t N_pop) const{ while(next_size < N_pop){ next_size += fronts[last_front++].size(); } - fronts.erase(fronts.begin() + last_front, fronts.end()); + fronts.erase(fronts.begin() + static_cast::difference_type>(last_front), fronts.end()); // Accept all members of first l-1 fronts for(size_t f=0; f nsga3::selection(population &R, size_t N_pop) const{ auto ext_points = find_extreme_points(R, fronts, translated_objectives); auto intercepts = find_intercepts(R, ext_points); auto norm_objs = normalize_objectives(translated_objectives, intercepts); - std::vector rps = generate_uniform_reference_points(nobj, divisions); + std::vector rps = generate_uniform_reference_points(nobj, m_divisions); associate_with_reference_points(rps, norm_objs, fronts); // Apply RP selection to final front until N_pop reached @@ -418,7 +416,7 @@ std::vector nsga3::selection(population &R, size_t N_pop) const{ rps[min_rp_idx].remove_candidate(selected_idx.value()); next.push_back(selected_idx.value()); }else{ - rps.erase(rps.begin() + min_rp_idx); + rps.erase(rps.begin() + static_cast::difference_type>(min_rp_idx)); } } @@ -428,7 +426,7 @@ std::vector nsga3::selection(population &R, size_t N_pop) const{ // Object serialization template void nsga3::serialize(Archive &ar, unsigned int) { - detail::archive(ar, ngen, cr, eta_c, mut, eta_mut, seed, m_verbosity, m_log); + detail::archive(ar, m_gen, m_cr, m_eta_c, m_mut, m_eta_mut, m_seed, m_verbosity, m_log); } } // namespace pagmo diff --git a/src/utils/multi_objective.cpp b/src/utils/multi_objective.cpp index 7adf3566d..2b7869a09 100644 --- a/src/utils/multi_objective.cpp +++ b/src/utils/multi_objective.cpp @@ -638,7 +638,16 @@ vector_double decompose_objectives(const vector_double &f, const vector_double & return {fd}; } -/* Gaussian Elimination */ +/// Gaussian Elimination +/** + * Solves the linear system Ax = b by performing Gaussian elimination + * on an augmented matrix A|b. + * + * @param + * @param + * @return + * @throws + */ vector_double gaussian_elimination(std::vector> A, const vector_double &b){ // N.B. Check dimensions and pivots if(A[0][0] == 0.0){ @@ -661,18 +670,32 @@ vector_double gaussian_elimination(std::vector> A, const vec } vector_double x(N); // Back substitution - for(int i=N-1; i>=0; i--){ + size_t i = N-1; + while(true){ for(size_t var=i+1; var::max(); @@ -686,7 +709,16 @@ double achievement(const vector_double &objs, const vector_double &weights){ return max_ratio; } -/* Distance from objective point to perpendicular intersection with reference point vector */ +/// Distance from objective point to perpendicular intersection with reference point vector */ +/** + * + * + * + * @param + * @param + * @return + * @throws + */ double perpendicular_distance(const std::vector &ref_point, const std::vector &obj_point){ double num = 0.0, denom = 0.0, sq_dist = 0.0; for(size_t i=0; i generate_reference_point_level( std::vector points; if(level == rp.dim()-1){ - rp[level] = 1.0*remain/total; + rp[level] = static_cast(remain)/static_cast(total); points.push_back(rp); }else{ for(size_t idx = 0; idx <= remain; idx++){ - rp[level] = 1.0*idx/total; + rp[level] = static_cast(idx)/static_cast(total); auto np = generate_reference_point_level(rp, remain - idx, level + 1, total); points.reserve(points.size() + np.size()); points.insert(points.end(), np.begin(), np.end()); @@ -67,7 +63,7 @@ void ReferencePoint::add_candidate(size_t index, double distance){ void ReferencePoint::remove_candidate(size_t index){ for(size_t idx=0; idx>::difference_type>(idx)); } } } @@ -142,7 +138,7 @@ std::optional ReferencePoint::random_candidate() const{ return choose_random_element>(candidates).first; } -size_t n_choose_k(unsigned n, unsigned k){ +size_t n_choose_k(size_t n, size_t k){ if(k == 0){ return 1u; } diff --git a/tests/nsga3.cpp b/tests/nsga3.cpp index 8926e6741..5e16e3ef4 100644 --- a/tests/nsga3.cpp +++ b/tests/nsga3.cpp @@ -13,7 +13,7 @@ #include #include #include -#include // gaussian_elimination +#include using namespace pagmo; @@ -69,23 +69,21 @@ BOOST_AUTO_TEST_CASE(nsga3_verify_uniform_reference_points){ } BOOST_AUTO_TEST_CASE(nsga3_test_translate_objectives){ - std::cout << "-==: nsga3_test_translate_objectives :==-\n"; + double tolerance = 1e-6; + std::vector t_first = {0.92084430016240049, 0.16973405319857038, 290.74413330194784}; + std::vector t_last = {1.7178358364136896, 109.71043974773266, 52.177816158337897}; + dtlz udp{1u, 10u, 3u}; population pop{udp, 52u, 23u}; nsga3 nsga3_alg{10u, 1.00, 30., 0.10, 20., 5u, 32u, false}; pop = nsga3_alg.evolve(pop); auto p0_obj = pop.get_f(); - std::cout << "-==: nsga3_test_translate_objectives pre translation :==-\n"; - for(size_t i=0; i < p0_obj.size(); i++){ - std::for_each(p0_obj[i].begin(), p0_obj[i].end(), [](const auto& elem){std::cout << elem << " "; }); - std::cout << std::endl; - } auto translated_objectives = nsga3_alg.translate_objectives(pop); - std::cout << "-==: nsga3_test_translate_objectives post translation :==-\n"; - for(size_t i=0; i < translated_objectives.size(); i++){ - std::for_each(translated_objectives[i].begin(), translated_objectives[i].end(), [](const auto& elem){std::cout << elem << " "; }); - std::cout << std::endl; + size_t t_size = translated_objectives.size(); + for(size_t i=0; i < translated_objectives[0].size(); i++){ + BOOST_CHECK_CLOSE(translated_objectives[0][i], t_first[i], tolerance); + BOOST_CHECK_CLOSE(translated_objectives[t_size-1][i], t_last[i], tolerance); } } @@ -94,16 +92,14 @@ BOOST_AUTO_TEST_CASE(nsga3_test_gaussian_elimination){ std::vector> A(3); std::vector b = {1.0, 1.0, 1.0}; - A[0] = {-1, 1, 2}; - A[1] = {2, 0, -3}; - A[2] = {5, 1, -2}; + A[0] = {-1, 1, 2}; + A[1] = { 2, 0, -3}; + A[2] = { 5, 1, -2}; std::vector x = gaussian_elimination(A, b); BOOST_CHECK_CLOSE(x[0], -0.4, 1e-8); BOOST_CHECK_CLOSE(x[1], 1.8, 1e-8); BOOST_CHECK_CLOSE(x[2], -0.6, 1e-8); - std::for_each(x.begin(), x.end(), [](const auto& i){std::cout << i << " ";}); - std::cout << std::endl; // Verify graceful error on ill-condition A[0][0] = 0.0; @@ -114,22 +110,30 @@ BOOST_AUTO_TEST_CASE(nsga3_test_find_extreme_points){ dtlz udp{1u, 10u, 3u}; population pop{udp, 52u, 23u}; nsga3 nsga3_alg{10u, 1.00, 30., 0.10, 20., 5u, 32u, false}; + std::vector ep_first = {228.71584572959793, 0.92448959747508574, 0.61400521336079161}; + std::vector ep_last = {0.092287013229137627, 0.0, 299.85225007963135}; + double tolerance = 1e-6; pop = nsga3_alg.evolve(pop); auto translated_objectives = nsga3_alg.translate_objectives(pop); auto fnds_res = fast_non_dominated_sorting(pop.get_f()); auto fronts = std::get<0>(fnds_res); auto ext_points = nsga3_alg.find_extreme_points(pop, fronts, translated_objectives); + size_t point_count = ext_points.size(); + size_t point_sz = ext_points[0].size(); - std::cout << "-==: extreme points :==-\n"; - //std::for_each(ext_points.begin(), ext_points.end(), [](const auto& elem){std::cout << elem << " "; }); - std::cout << std::endl; + for(size_t idx=0; idx < point_sz; idx++){ + BOOST_CHECK_CLOSE(ext_points[0][idx], ep_first[idx], tolerance); + BOOST_CHECK_CLOSE(ext_points[point_count-1][idx], ep_last[idx], tolerance); + } } BOOST_AUTO_TEST_CASE(nsga3_test_find_intercepts){ dtlz udp{1u, 10u, 3u}; population pop{udp, 52u, 23u}; nsga3 nsga3_alg{10u, 1.00, 30., 0.10, 20., 5u, 32u, false}; + std::vector intercept_values = {230.13712800033696, 223.90511605342394, 299.97254170821623}; + double tolerance = 1e-6; pop = nsga3_alg.evolve(pop); auto translated_objectives = nsga3_alg.translate_objectives(pop); @@ -138,15 +142,18 @@ BOOST_AUTO_TEST_CASE(nsga3_test_find_intercepts){ auto ext_points = nsga3_alg.find_extreme_points(pop, fronts, translated_objectives); auto intercepts = nsga3_alg.find_intercepts(pop, ext_points); - std::cout << "-==: intercepts :==-\n"; - std::for_each(intercepts.begin(), intercepts.end(), [](const auto& elem){std::cout << elem << " "; }); - std::cout << std::endl; + for(size_t idx=0; idx < intercepts.size(); idx++){ + BOOST_CHECK_CLOSE(intercepts[idx], intercept_values[idx], tolerance); + } } BOOST_AUTO_TEST_CASE(nsga3_test_normalize_objectives){ dtlz udp{1u, 10u, 3u}; population pop{udp, 52u, 23u}; nsga3 nsga3_alg{10u, 1.00, 30., 0.10, 20., 5u, 32u, false}; + std::vector norm_first = {0.0040012852691941663, 0.00075806241585865187, 0.96923582287326526}; + std::vector norm_last = {0.0074644011218006267, 0.48998630170449375, 0.173941974359411}; + double tolerance = 1e-6; pop = nsga3_alg.evolve(pop); auto translated_objectives = nsga3_alg.translate_objectives(pop); @@ -155,16 +162,14 @@ BOOST_AUTO_TEST_CASE(nsga3_test_normalize_objectives){ auto ext_points = nsga3_alg.find_extreme_points(pop, fronts, translated_objectives); auto intercepts = nsga3_alg.find_intercepts(pop, ext_points); auto norm_objs = nsga3_alg.normalize_objectives(translated_objectives, intercepts); - for(const auto &obj_f: norm_objs){ - std::for_each(obj_f.begin(), obj_f.end(), [](const auto& elem){std::cout << elem << " "; }); - std::cout << std::endl; + size_t obj_count = norm_objs.size(); + size_t obj_len = norm_objs[0].size(); + for(size_t idx=0; idx < obj_len; idx++){ + BOOST_CHECK_CLOSE(norm_objs[0][idx], norm_first[idx], tolerance); + BOOST_CHECK_CLOSE(norm_objs[obj_count-1][idx], norm_last[idx], tolerance); } } -//BOOST_AUTO_TEST_CASE(nsga3_test_associate_reference_points){ -// nsga3 n{10u, 0.95, 10., 0.01, 50., 32u, false}; -//} - BOOST_AUTO_TEST_CASE(nsga3_serialization_test){ double close_distance = 1e-8; problem prob{zdt{1u, 30u}};