From 847d2dfdb97509fdcff2e4bb8130e384a590a8e5 Mon Sep 17 00:00:00 2001 From: Evgenii Vikulov Date: Mon, 5 Jun 2023 14:43:53 +0900 Subject: [PATCH 1/5] Added fitness batch evaluation to enable multi threading --- CMakeLists.txt | 2 +- include/pagmo/algorithms/de.hpp | 17 ++ src/algorithms/de.cpp | 379 +++++++++++++++++++------------- tests/de.cpp | 33 +++ 4 files changed, 279 insertions(+), 152 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 2364363b4..e546841bf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -22,7 +22,7 @@ message(STATUS "System name: ${CMAKE_SYSTEM_NAME}") enable_testing() # Build option: enable the test suite. -option(PAGMO_BUILD_TESTS "Build the test suite." OFF) +option(PAGMO_BUILD_TESTS "Build the test suite." ON) # Build option: enable the benchmarks. option(PAGMO_BUILD_BENCHMARKS "Build the benchmarks." OFF) diff --git a/include/pagmo/algorithms/de.hpp b/include/pagmo/algorithms/de.hpp index 29fc02e39..8715fd5d8 100644 --- a/include/pagmo/algorithms/de.hpp +++ b/include/pagmo/algorithms/de.hpp @@ -33,7 +33,10 @@ see https://www.gnu.org/licenses/. */ #include #include +#include + #include +#include #include #include #include @@ -159,6 +162,10 @@ class PAGMO_DLL_PUBLIC de { return m_gen; } + + // Sets the bfe + void set_bfe(const bfe &b); + /// Algorithm name /** * One of the optional methods of any user-defined algorithm (UDA). @@ -184,6 +191,15 @@ class PAGMO_DLL_PUBLIC de } private: + vector_double mutate(const population &pop, population::size_type i, const vector_double &gbIter, + std::uniform_real_distribution drng, + std::uniform_int_distribution c_idx) const; + + void update_pop(population &pop, const vector_double &newfitness, population::size_type i, + std::vector &fit, vector_double &gbfit, vector_double &gbX, + std::vector &popnew, const std::vector &popold, + const vector_double &tmp) const; + // Object serialization friend class boost::serialization::access; template @@ -199,6 +215,7 @@ class PAGMO_DLL_PUBLIC de unsigned m_seed; unsigned m_verbosity; mutable log_type m_log; + boost::optional m_bfe; }; } // namespace pagmo diff --git a/src/algorithms/de.cpp b/src/algorithms/de.cpp index bdbb819ad..38db3c6f0 100644 --- a/src/algorithms/de.cpp +++ b/src/algorithms/de.cpp @@ -48,6 +48,174 @@ see https://www.gnu.org/licenses/. */ namespace pagmo { +vector_double de::mutate(const population &pop, population::size_type i, const vector_double &gbIter, + std::uniform_real_distribution drng, + std::uniform_int_distribution c_idx) const +{ + const auto &prob = pop.get_problem(); // This is a const reference, so using set_seed for example will not be + // allowed + auto dim = prob.get_nx(); + auto NP = pop.size(); + auto popold = pop.get_x(); + std::vector r(5); // indexes of 5 selected population members + vector_double tmp(dim); + + /*-----We select at random 5 indexes from the population---------------------------------*/ + std::vector idxs(NP); + std::iota(idxs.begin(), idxs.end(), vector_double::size_type(0u)); + for (auto j = 0u; j < 5u; ++j) { // Durstenfeld's algorithm to select 5 indexes at random + auto idx = std::uniform_int_distribution(0u, NP - 1u - j)(m_e); + r[j] = idxs[idx]; + std::swap(idxs[idx], idxs[NP - 1u - j]); + } + + /*-------DE/best/1/exp--------------------------------------------------------------------*/ + /*-------The oldest DE variant but still not bad. However, we have found several---------*/ + /*-------optimization problems where misconvergence occurs.-------------------------------*/ + if (m_variant == 1u) { + tmp = popold[i]; + auto n = c_idx(m_e); + auto L = 0u; + do { + tmp[n] = gbIter[n] + m_F * (popold[r[1]][n] - popold[r[2]][n]); + n = (n + 1u) % dim; + ++L; + } while ((drng(m_e) < m_CR) && (L < dim)); + } + + /*-------DE/rand/1/exp-------------------------------------------------------------------*/ + /*-------This is one of my favourite strategies. It works especially well when the-------*/ + /*-------"gbIter[]"-schemes experience misconvergence. Try e.g. m_F=0.7 and m_CR=0.5---------*/ + /*-------as a first guess.---------------------------------------------------------------*/ + else if (m_variant == 2u) { + tmp = popold[i]; + auto n = c_idx(m_e); + decltype(dim) L = 0u; + do { + tmp[n] = popold[r[0]][n] + m_F * (popold[r[1]][n] - popold[r[2]][n]); + n = (n + 1u) % dim; + ++L; + } while ((drng(m_e) < m_CR) && (L < dim)); + } + /*-------DE/rand-to-best/1/exp-----------------------------------------------------------*/ + /*-------This variant seems to be one of the best strategies. Try m_F=0.85 and m_CR=1.------*/ + /*-------If you get misconvergence try to increase NP. If this doesn't help you----------*/ + /*-------should play around with all three control variables.----------------------------*/ + else if (m_variant == 3u) { + tmp = popold[i]; + auto n = c_idx(m_e); + auto L = 0u; + do { + tmp[n] = tmp[n] + m_F * (gbIter[n] - tmp[n]) + m_F * (popold[r[0]][n] - popold[r[1]][n]); + n = (n + 1u) % dim; + ++L; + } while ((drng(m_e) < m_CR) && (L < dim)); + } + /*-------DE/best/2/exp is another powerful variant worth trying--------------------------*/ + else if (m_variant == 4u) { + tmp = popold[i]; + auto n = c_idx(m_e); + auto L = 0u; + do { + tmp[n] = gbIter[n] + (popold[r[0]][n] + popold[r[1]][n] - popold[r[2]][n] - popold[r[3]][n]) * m_F; + n = (n + 1u) % dim; + ++L; + } while ((drng(m_e) < m_CR) && (L < dim)); + } + /*-------DE/rand/2/exp seems to be a robust optimizer for many functions-------------------*/ + else if (m_variant == 5u) { + tmp = popold[i]; + auto n = c_idx(m_e); + auto L = 0u; + do { + tmp[n] = popold[r[4]][n] + (popold[r[0]][n] + popold[r[1]][n] - popold[r[2]][n] - popold[r[3]][n]) * m_F; + n = (n + 1u) % dim; + ++L; + } while ((drng(m_e) < m_CR) && (L < dim)); + } + + /*=======Essentially same strategies but BINOMIAL CROSSOVER===============================*/ + /*-------DE/best/1/bin--------------------------------------------------------------------*/ + else if (m_variant == 6u) { + tmp = popold[i]; + auto n = c_idx(m_e); + for (decltype(dim) L = 0u; L < dim; ++L) { /* perform Dc binomial trials */ + if ((drng(m_e) < m_CR) || L + 1u == dim) { /* change at least one parameter */ + tmp[n] = gbIter[n] + m_F * (popold[r[1]][n] - popold[r[2]][n]); + } + n = (n + 1u) % dim; + } + } + /*-------DE/rand/1/bin-------------------------------------------------------------------*/ + else if (m_variant == 7u) { + tmp = popold[i]; + auto n = c_idx(m_e); + for (decltype(dim) L = 0u; L < dim; ++L) { /* perform Dc binomial trials */ + if ((drng(m_e) < m_CR) || L + 1u == dim) { /* change at least one parameter */ + tmp[n] = popold[r[0]][n] + m_F * (popold[r[1]][n] - popold[r[2]][n]); + } + n = (n + 1u) % dim; + } + } + /*-------DE/rand-to-best/1/bin-----------------------------------------------------------*/ + else if (m_variant == 8u) { + tmp = popold[i]; + auto n = c_idx(m_e); + for (decltype(dim) L = 0u; L < dim; ++L) { /* perform Dc binomial trials */ + if ((drng(m_e) < m_CR) || L + 1u == dim) { /* change at least one parameter */ + tmp[n] = tmp[n] + m_F * (gbIter[n] - tmp[n]) + m_F * (popold[r[0]][n] - popold[r[1]][n]); + } + n = (n + 1u) % dim; + } + } + /*-------DE/best/2/bin--------------------------------------------------------------------*/ + else if (m_variant == 9u) { + tmp = popold[i]; + auto n = c_idx(m_e); + for (decltype(dim) L = 0u; L < dim; ++L) { /* perform Dc binomial trials */ + if ((drng(m_e) < m_CR) || L + 1u == dim) { /* change at least one parameter */ + tmp[n] = gbIter[n] + (popold[r[0]][n] + popold[r[1]][n] - popold[r[2]][n] - popold[r[3]][n]) * m_F; + } + n = (n + 1u) % dim; + } + } + /*-------DE/rand/2/bin--------------------------------------------------------------------*/ + else if (m_variant == 10u) { + tmp = popold[i]; + auto n = c_idx(m_e); + for (decltype(dim) L = 0u; L < dim; ++L) { /* perform Dc binomial trials */ + if ((drng(m_e) < m_CR) || L + 1u == dim) { /* change at least one parameter */ + tmp[n] + = popold[r[4]][n] + (popold[r[0]][n] + popold[r[1]][n] - popold[r[2]][n] - popold[r[3]][n]) * m_F; + } + n = (n + 1u) % dim; + } + } + + return tmp; +} + +void de::update_pop(population &pop, const vector_double &newfitness, population::size_type i, + std::vector &fit, vector_double &gbfit, vector_double &gbX, + std::vector &popnew, const std::vector &popold, + const vector_double &tmp) const +{ + if (newfitness[0] <= fit[i][0]) { /* improved objective function value ? */ + fit[i] = newfitness; + popnew[i] = tmp; + // updates the individual in pop (avoiding to recompute the objective function) + pop.set_xf(i, popnew[i], newfitness); + + if (newfitness[0] <= gbfit[0]) { + /* if so...*/ + gbfit = newfitness; /* reset gbfit to new low...*/ + gbX = popnew[i]; + } + } else { + popnew[i] = popold[i]; + } +} + de::de(unsigned gen, double F, double CR, unsigned variant, double ftol, double xtol, unsigned seed) : m_gen(gen), m_F(F), m_CR(CR), m_variant(variant), m_Ftol(ftol), m_xtol(xtol), m_e(seed), m_seed(seed), m_verbosity(0u), m_log() @@ -115,8 +283,6 @@ population de::evolve(population pop) const // No throws, all valid: we clear the logs m_log.clear(); - // Some vectors used during evolution are declared. - vector_double tmp(dim); // contains the mutated candidate std::uniform_real_distribution drng(0., 1.); // to generate a number in [0, 1) std::uniform_int_distribution c_idx( 0u, dim - 1u); // to generate a random index for the chromosome @@ -133,167 +299,69 @@ population de::evolve(population pop) const auto gbfit = fit[best_idx]; // the best decision vector of a generation auto gbIter = gbX; - std::vector r(5); // indexes of 5 selected population members + + // Main DE iterations for (decltype(m_gen) gen = 1u; gen <= m_gen; ++gen) { - // Start of the loop through the population - for (decltype(NP) i = 0u; i < NP; ++i) { - /*-----We select at random 5 indexes from the population---------------------------------*/ - std::vector idxs(NP); - std::iota(idxs.begin(), idxs.end(), vector_double::size_type(0u)); - for (auto j = 0u; j < 5u; ++j) { // Durstenfeld's algorithm to select 5 indexes at random - auto idx = std::uniform_int_distribution(0u, NP - 1u - j)(m_e); - r[j] = idxs[idx]; - std::swap(idxs[idx], idxs[NP - 1u - j]); - } + + if (m_bfe) { + // bfe is available: + vector_double trial(NP * dim); + decltype(trial.size()) pos = 0u; - /*-------DE/best/1/exp--------------------------------------------------------------------*/ - /*-------The oldest DE variant but still not bad. However, we have found several---------*/ - /*-------optimization problems where misconvergence occurs.-------------------------------*/ - if (m_variant == 1u) { - tmp = popold[i]; - auto n = c_idx(m_e); - auto L = 0u; - do { - tmp[n] = gbIter[n] + m_F * (popold[r[1]][n] - popold[r[2]][n]); - n = (n + 1u) % dim; - ++L; - } while ((drng(m_e) < m_CR) && (L < dim)); - } + // Start of the loop through the population + for (decltype(NP) i = 0u; i < NP; ++i) { - /*-------DE/rand/1/exp-------------------------------------------------------------------*/ - /*-------This is one of my favourite strategies. It works especially well when the-------*/ - /*-------"gbIter[]"-schemes experience misconvergence. Try e.g. m_F=0.7 and m_CR=0.5---------*/ - /*-------as a first guess.---------------------------------------------------------------*/ - else if (m_variant == 2u) { - tmp = popold[i]; - auto n = c_idx(m_e); - decltype(dim) L = 0u; - do { - tmp[n] = popold[r[0]][n] + m_F * (popold[r[1]][n] - popold[r[2]][n]); - n = (n + 1u) % dim; - ++L; - } while ((drng(m_e) < m_CR) && (L < dim)); - } - /*-------DE/rand-to-best/1/exp-----------------------------------------------------------*/ - /*-------This variant seems to be one of the best strategies. Try m_F=0.85 and m_CR=1.------*/ - /*-------If you get misconvergence try to increase NP. If this doesn't help you----------*/ - /*-------should play around with all three control variables.----------------------------*/ - else if (m_variant == 3u) { - tmp = popold[i]; - auto n = c_idx(m_e); - auto L = 0u; - do { - tmp[n] = tmp[n] + m_F * (gbIter[n] - tmp[n]) + m_F * (popold[r[0]][n] - popold[r[1]][n]); - n = (n + 1u) % dim; - ++L; - } while ((drng(m_e) < m_CR) && (L < dim)); - } - /*-------DE/best/2/exp is another powerful variant worth trying--------------------------*/ - else if (m_variant == 4u) { - tmp = popold[i]; - auto n = c_idx(m_e); - auto L = 0u; - do { - tmp[n] = gbIter[n] + (popold[r[0]][n] + popold[r[1]][n] - popold[r[2]][n] - popold[r[3]][n]) * m_F; - n = (n + 1u) % dim; - ++L; - } while ((drng(m_e) < m_CR) && (L < dim)); - } - /*-------DE/rand/2/exp seems to be a robust optimizer for many functions-------------------*/ - else if (m_variant == 5u) { - tmp = popold[i]; - auto n = c_idx(m_e); - auto L = 0u; - do { - tmp[n] = popold[r[4]][n] - + (popold[r[0]][n] + popold[r[1]][n] - popold[r[2]][n] - popold[r[3]][n]) * m_F; - n = (n + 1u) % dim; - ++L; - } while ((drng(m_e) < m_CR) && (L < dim)); - } + auto tmp = mutate(pop, i, gbIter, drng, c_idx); - /*=======Essentially same strategies but BINOMIAL CROSSOVER===============================*/ - /*-------DE/best/1/bin--------------------------------------------------------------------*/ - else if (m_variant == 6u) { - tmp = popold[i]; - auto n = c_idx(m_e); - for (decltype(dim) L = 0u; L < dim; ++L) { /* perform Dc binomial trials */ - if ((drng(m_e) < m_CR) || L + 1u == dim) { /* change at least one parameter */ - tmp[n] = gbIter[n] + m_F * (popold[r[1]][n] - popold[r[2]][n]); - } - n = (n + 1u) % dim; - } - } - /*-------DE/rand/1/bin-------------------------------------------------------------------*/ - else if (m_variant == 7u) { - tmp = popold[i]; - auto n = c_idx(m_e); - for (decltype(dim) L = 0u; L < dim; ++L) { /* perform Dc binomial trials */ - if ((drng(m_e) < m_CR) || L + 1u == dim) { /* change at least one parameter */ - tmp[n] = popold[r[0]][n] + m_F * (popold[r[1]][n] - popold[r[2]][n]); - } - n = (n + 1u) % dim; - } - } - /*-------DE/rand-to-best/1/bin-----------------------------------------------------------*/ - else if (m_variant == 8u) { - tmp = popold[i]; - auto n = c_idx(m_e); - for (decltype(dim) L = 0u; L < dim; ++L) { /* perform Dc binomial trials */ - if ((drng(m_e) < m_CR) || L + 1u == dim) { /* change at least one parameter */ - tmp[n] = tmp[n] + m_F * (gbIter[n] - tmp[n]) + m_F * (popold[r[0]][n] - popold[r[1]][n]); - } - n = (n + 1u) % dim; - } - } - /*-------DE/best/2/bin--------------------------------------------------------------------*/ - else if (m_variant == 9u) { - tmp = popold[i]; - auto n = c_idx(m_e); - for (decltype(dim) L = 0u; L < dim; ++L) { /* perform Dc binomial trials */ - if ((drng(m_e) < m_CR) || L + 1u == dim) { /* change at least one parameter */ - tmp[n] - = gbIter[n] + (popold[r[0]][n] + popold[r[1]][n] - popold[r[2]][n] - popold[r[3]][n]) * m_F; - } - n = (n + 1u) % dim; + // Trial mutation now in tmp. force feasibility and see how good this choice really was. + // a) feasibility + // detail::force_bounds_reflection(tmp, lb, ub); // TODO: check if this choice is better + detail::force_bounds_random(tmp, lb, ub, m_e); + + for (decltype(tmp.size()) ii = 0u; ii < tmp.size(); ++ii) { + trial[pos] = tmp[ii]; + ++pos; } } - /*-------DE/rand/2/bin--------------------------------------------------------------------*/ - else if (m_variant == 10u) { - tmp = popold[i]; - auto n = c_idx(m_e); - for (decltype(dim) L = 0u; L < dim; ++L) { /* perform Dc binomial trials */ - if ((drng(m_e) < m_CR) || L + 1u == dim) { /* change at least one parameter */ - tmp[n] = popold[r[4]][n] - + (popold[r[0]][n] + popold[r[1]][n] - popold[r[2]][n] - popold[r[3]][n]) * m_F; - } - n = (n + 1u) % dim; + + vector_double tmp(dim); + vector_double newfitness(prob_f_dimension); + auto fitnesses = (*m_bfe)(prob, trial); + decltype(tmp.size()) pos_dim = 0u; + decltype(newfitness.size()) pos_fit = 0u; + + for (decltype(NP) i = 0u; i < NP; ++i) { + for (decltype(dim) ii_dim = 0u; ii_dim < dim; ++ii_dim) { + tmp[ii_dim] = trial[pos_dim]; + ++pos_dim; } - } - // Trial mutation now in tmp. force feasibility and see how good this choice really was. - // a) feasibility - // detail::force_bounds_reflection(tmp, lb, ub); // TODO: check if this choice is better - detail::force_bounds_random(tmp, lb, ub, m_e); - // b) how good? - auto newfitness = prob.fitness(tmp); /* Evaluates tmp[] */ - if (newfitness[0] <= fit[i][0]) { /* improved objective function value ? */ - fit[i] = newfitness; - popnew[i] = tmp; - // updates the individual in pop (avoiding to recompute the objective function) - pop.set_xf(i, popnew[i], newfitness); - - if (newfitness[0] <= gbfit[0]) { - /* if so...*/ - gbfit = newfitness; /* reset gbfit to new low...*/ - gbX = popnew[i]; + for (decltype(prob_f_dimension) ii_f = 0u; ii_f < prob_f_dimension; ++ii_f) { + newfitness[ii_f] = fitnesses[pos_fit]; + ++pos_fit; } - } else { - popnew[i] = popold[i]; + + update_pop(pop, newfitness, i, fit, gbfit, gbX, popnew, popold, tmp); } - } // End of one generation + + } else { + + for (decltype(NP) i = 0u; i < NP; ++i) { + + auto tmp = mutate(pop, i, gbIter, drng, c_idx); + // Trial mutation now in tmp. force feasibility and see how good this choice really was. + // a) feasibility + // detail::force_bounds_reflection(tmp, lb, ub); // TODO: check if this choice is better + detail::force_bounds_random(tmp, lb, ub, m_e); + // b) how good? + auto newfitness = prob.fitness(tmp); /* Evaluates tmp[] */ + + update_pop(pop, newfitness, i, fit, gbfit, gbX, popnew, popold, tmp); + } // End of one generation + } + /* Save best population member of current iteration */ gbIter = gbX; /* swap population arrays. New generation becomes old one */ @@ -363,6 +431,15 @@ void de::set_seed(unsigned seed) m_seed = seed; } +/// Sets the batch function evaluation scheme +/** + * @param b batch function evaluation object + */ +void de::set_bfe(const bfe &b) +{ + m_bfe = b; +} + /// Extra info /** * One of the optional methods of any user-defined algorithm (UDA). diff --git a/tests/de.cpp b/tests/de.cpp index c79083202..8ddeb566c 100644 --- a/tests/de.cpp +++ b/tests/de.cpp @@ -32,6 +32,7 @@ see https://www.gnu.org/licenses/. */ #include #include +#include #include #include @@ -174,3 +175,35 @@ BOOST_AUTO_TEST_CASE(de_serialization_test) BOOST_CHECK_CLOSE(std::get<4>(before_log[i]), std::get<4>(after_log[i]), 1e-8); } } + + +BOOST_AUTO_TEST_CASE(bfe_usage_test) +{ + + std::chrono::steady_clock::time_point begin_1 = std::chrono::steady_clock::now(); + population pop{rosenbrock{10u}, 40u, 23u}; + de user_algo{1000000u, 0.7, 0.5, 2u, 1e-3, 1e-50, 23u}; + user_algo.set_verbosity(1u); + user_algo.set_seed(23u); + user_algo.set_bfe(bfe{}); + pop = user_algo.evolve(pop); + std::chrono::steady_clock::time_point end_1 = std::chrono::steady_clock::now(); + BOOST_CHECK(user_algo.get_log().size() < 6000u); + + std::chrono::steady_clock::time_point begin_2 = std::chrono::steady_clock::now(); + population pop_2{rosenbrock{10u}, 40u, 23u}; + de user_algo_2{1000000u, 0.7, 0.5, 2u, 1e-3, 1e-50, 23u}; + user_algo_2.set_verbosity(1u); + pop_2 = user_algo_2.evolve(pop_2); + std::chrono::steady_clock::time_point end_2 = std::chrono::steady_clock::now(); + BOOST_CHECK(user_algo_2.get_log().size() < 6000u); + + BOOST_CHECK_EQUAL(pop.champion_f()[0], pop_2.champion_f()[0]); + + + std::cout << "Time difference = " << std::chrono::duration_cast(end_1 - begin_1).count() + << "[ms]" << std::endl; + + std::cout << "Time difference = " << std::chrono::duration_cast(end_2 - begin_2).count() + << "[ms]" << std::endl; +} From 9859b9f81b81281f68dc0cc43c442c2e7194f682 Mon Sep 17 00:00:00 2001 From: Evgenii Vikulov Date: Mon, 5 Jun 2023 18:34:59 +0900 Subject: [PATCH 2/5] Removed expensive copy --- include/pagmo/algorithms/de.hpp | 3 ++- src/algorithms/de.cpp | 8 ++++---- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/include/pagmo/algorithms/de.hpp b/include/pagmo/algorithms/de.hpp index 8715fd5d8..86614b7e1 100644 --- a/include/pagmo/algorithms/de.hpp +++ b/include/pagmo/algorithms/de.hpp @@ -193,7 +193,8 @@ class PAGMO_DLL_PUBLIC de private: vector_double mutate(const population &pop, population::size_type i, const vector_double &gbIter, std::uniform_real_distribution drng, - std::uniform_int_distribution c_idx) const; + std::uniform_int_distribution c_idx, + const std::vector& popold) const; void update_pop(population &pop, const vector_double &newfitness, population::size_type i, std::vector &fit, vector_double &gbfit, vector_double &gbX, diff --git a/src/algorithms/de.cpp b/src/algorithms/de.cpp index 38db3c6f0..8c15fe249 100644 --- a/src/algorithms/de.cpp +++ b/src/algorithms/de.cpp @@ -50,13 +50,13 @@ namespace pagmo vector_double de::mutate(const population &pop, population::size_type i, const vector_double &gbIter, std::uniform_real_distribution drng, - std::uniform_int_distribution c_idx) const + std::uniform_int_distribution c_idx, + const std::vector& popold) const { const auto &prob = pop.get_problem(); // This is a const reference, so using set_seed for example will not be // allowed auto dim = prob.get_nx(); auto NP = pop.size(); - auto popold = pop.get_x(); std::vector r(5); // indexes of 5 selected population members vector_double tmp(dim); @@ -313,7 +313,7 @@ population de::evolve(population pop) const // Start of the loop through the population for (decltype(NP) i = 0u; i < NP; ++i) { - auto tmp = mutate(pop, i, gbIter, drng, c_idx); + auto tmp = mutate(pop, i, gbIter, drng, c_idx, popold); // Trial mutation now in tmp. force feasibility and see how good this choice really was. // a) feasibility @@ -350,7 +350,7 @@ population de::evolve(population pop) const for (decltype(NP) i = 0u; i < NP; ++i) { - auto tmp = mutate(pop, i, gbIter, drng, c_idx); + auto tmp = mutate(pop, i, gbIter, drng, c_idx, popold); // Trial mutation now in tmp. force feasibility and see how good this choice really was. // a) feasibility // detail::force_bounds_reflection(tmp, lb, ub); // TODO: check if this choice is better From 3c2ba13523f0097bd478c9de1e5af16fd83618e0 Mon Sep 17 00:00:00 2001 From: Evgenii Vikulov Date: Tue, 6 Jun 2023 11:57:37 +0900 Subject: [PATCH 3/5] tidy up --- src/algorithms/de.cpp | 345 +++++++++++++++++++++--------------------- tests/de.cpp | 24 +-- 2 files changed, 177 insertions(+), 192 deletions(-) diff --git a/src/algorithms/de.cpp b/src/algorithms/de.cpp index 8c15fe249..9ec7e6c1a 100644 --- a/src/algorithms/de.cpp +++ b/src/algorithms/de.cpp @@ -45,177 +45,12 @@ see https://www.gnu.org/licenses/. */ #include #include -namespace pagmo -{ - -vector_double de::mutate(const population &pop, population::size_type i, const vector_double &gbIter, - std::uniform_real_distribution drng, - std::uniform_int_distribution c_idx, - const std::vector& popold) const -{ - const auto &prob = pop.get_problem(); // This is a const reference, so using set_seed for example will not be - // allowed - auto dim = prob.get_nx(); - auto NP = pop.size(); - std::vector r(5); // indexes of 5 selected population members - vector_double tmp(dim); +// NOTE: apparently this must be included *after* +// the other serialization headers. +#include - /*-----We select at random 5 indexes from the population---------------------------------*/ - std::vector idxs(NP); - std::iota(idxs.begin(), idxs.end(), vector_double::size_type(0u)); - for (auto j = 0u; j < 5u; ++j) { // Durstenfeld's algorithm to select 5 indexes at random - auto idx = std::uniform_int_distribution(0u, NP - 1u - j)(m_e); - r[j] = idxs[idx]; - std::swap(idxs[idx], idxs[NP - 1u - j]); - } - - /*-------DE/best/1/exp--------------------------------------------------------------------*/ - /*-------The oldest DE variant but still not bad. However, we have found several---------*/ - /*-------optimization problems where misconvergence occurs.-------------------------------*/ - if (m_variant == 1u) { - tmp = popold[i]; - auto n = c_idx(m_e); - auto L = 0u; - do { - tmp[n] = gbIter[n] + m_F * (popold[r[1]][n] - popold[r[2]][n]); - n = (n + 1u) % dim; - ++L; - } while ((drng(m_e) < m_CR) && (L < dim)); - } - - /*-------DE/rand/1/exp-------------------------------------------------------------------*/ - /*-------This is one of my favourite strategies. It works especially well when the-------*/ - /*-------"gbIter[]"-schemes experience misconvergence. Try e.g. m_F=0.7 and m_CR=0.5---------*/ - /*-------as a first guess.---------------------------------------------------------------*/ - else if (m_variant == 2u) { - tmp = popold[i]; - auto n = c_idx(m_e); - decltype(dim) L = 0u; - do { - tmp[n] = popold[r[0]][n] + m_F * (popold[r[1]][n] - popold[r[2]][n]); - n = (n + 1u) % dim; - ++L; - } while ((drng(m_e) < m_CR) && (L < dim)); - } - /*-------DE/rand-to-best/1/exp-----------------------------------------------------------*/ - /*-------This variant seems to be one of the best strategies. Try m_F=0.85 and m_CR=1.------*/ - /*-------If you get misconvergence try to increase NP. If this doesn't help you----------*/ - /*-------should play around with all three control variables.----------------------------*/ - else if (m_variant == 3u) { - tmp = popold[i]; - auto n = c_idx(m_e); - auto L = 0u; - do { - tmp[n] = tmp[n] + m_F * (gbIter[n] - tmp[n]) + m_F * (popold[r[0]][n] - popold[r[1]][n]); - n = (n + 1u) % dim; - ++L; - } while ((drng(m_e) < m_CR) && (L < dim)); - } - /*-------DE/best/2/exp is another powerful variant worth trying--------------------------*/ - else if (m_variant == 4u) { - tmp = popold[i]; - auto n = c_idx(m_e); - auto L = 0u; - do { - tmp[n] = gbIter[n] + (popold[r[0]][n] + popold[r[1]][n] - popold[r[2]][n] - popold[r[3]][n]) * m_F; - n = (n + 1u) % dim; - ++L; - } while ((drng(m_e) < m_CR) && (L < dim)); - } - /*-------DE/rand/2/exp seems to be a robust optimizer for many functions-------------------*/ - else if (m_variant == 5u) { - tmp = popold[i]; - auto n = c_idx(m_e); - auto L = 0u; - do { - tmp[n] = popold[r[4]][n] + (popold[r[0]][n] + popold[r[1]][n] - popold[r[2]][n] - popold[r[3]][n]) * m_F; - n = (n + 1u) % dim; - ++L; - } while ((drng(m_e) < m_CR) && (L < dim)); - } - - /*=======Essentially same strategies but BINOMIAL CROSSOVER===============================*/ - /*-------DE/best/1/bin--------------------------------------------------------------------*/ - else if (m_variant == 6u) { - tmp = popold[i]; - auto n = c_idx(m_e); - for (decltype(dim) L = 0u; L < dim; ++L) { /* perform Dc binomial trials */ - if ((drng(m_e) < m_CR) || L + 1u == dim) { /* change at least one parameter */ - tmp[n] = gbIter[n] + m_F * (popold[r[1]][n] - popold[r[2]][n]); - } - n = (n + 1u) % dim; - } - } - /*-------DE/rand/1/bin-------------------------------------------------------------------*/ - else if (m_variant == 7u) { - tmp = popold[i]; - auto n = c_idx(m_e); - for (decltype(dim) L = 0u; L < dim; ++L) { /* perform Dc binomial trials */ - if ((drng(m_e) < m_CR) || L + 1u == dim) { /* change at least one parameter */ - tmp[n] = popold[r[0]][n] + m_F * (popold[r[1]][n] - popold[r[2]][n]); - } - n = (n + 1u) % dim; - } - } - /*-------DE/rand-to-best/1/bin-----------------------------------------------------------*/ - else if (m_variant == 8u) { - tmp = popold[i]; - auto n = c_idx(m_e); - for (decltype(dim) L = 0u; L < dim; ++L) { /* perform Dc binomial trials */ - if ((drng(m_e) < m_CR) || L + 1u == dim) { /* change at least one parameter */ - tmp[n] = tmp[n] + m_F * (gbIter[n] - tmp[n]) + m_F * (popold[r[0]][n] - popold[r[1]][n]); - } - n = (n + 1u) % dim; - } - } - /*-------DE/best/2/bin--------------------------------------------------------------------*/ - else if (m_variant == 9u) { - tmp = popold[i]; - auto n = c_idx(m_e); - for (decltype(dim) L = 0u; L < dim; ++L) { /* perform Dc binomial trials */ - if ((drng(m_e) < m_CR) || L + 1u == dim) { /* change at least one parameter */ - tmp[n] = gbIter[n] + (popold[r[0]][n] + popold[r[1]][n] - popold[r[2]][n] - popold[r[3]][n]) * m_F; - } - n = (n + 1u) % dim; - } - } - /*-------DE/rand/2/bin--------------------------------------------------------------------*/ - else if (m_variant == 10u) { - tmp = popold[i]; - auto n = c_idx(m_e); - for (decltype(dim) L = 0u; L < dim; ++L) { /* perform Dc binomial trials */ - if ((drng(m_e) < m_CR) || L + 1u == dim) { /* change at least one parameter */ - tmp[n] - = popold[r[4]][n] + (popold[r[0]][n] + popold[r[1]][n] - popold[r[2]][n] - popold[r[3]][n]) * m_F; - } - n = (n + 1u) % dim; - } - } - - return tmp; -} - -void de::update_pop(population &pop, const vector_double &newfitness, population::size_type i, - std::vector &fit, vector_double &gbfit, vector_double &gbX, - std::vector &popnew, const std::vector &popold, - const vector_double &tmp) const +namespace pagmo { - if (newfitness[0] <= fit[i][0]) { /* improved objective function value ? */ - fit[i] = newfitness; - popnew[i] = tmp; - // updates the individual in pop (avoiding to recompute the objective function) - pop.set_xf(i, popnew[i], newfitness); - - if (newfitness[0] <= gbfit[0]) { - /* if so...*/ - gbfit = newfitness; /* reset gbfit to new low...*/ - gbX = popnew[i]; - } - } else { - popnew[i] = popold[i]; - } -} - de::de(unsigned gen, double F, double CR, unsigned variant, double ftol, double xtol, unsigned seed) : m_gen(gen), m_F(F), m_CR(CR), m_variant(variant), m_Ftol(ftol), m_xtol(xtol), m_e(seed), m_seed(seed), m_verbosity(0u), m_log() @@ -300,8 +135,6 @@ population de::evolve(population pop) const // the best decision vector of a generation auto gbIter = gbX; - - // Main DE iterations for (decltype(m_gen) gen = 1u; gen <= m_gen; ++gen) { @@ -421,6 +254,174 @@ population de::evolve(population pop) const return pop; } +vector_double de::mutate(const population &pop, population::size_type i, const vector_double &gbIter, + std::uniform_real_distribution drng, + std::uniform_int_distribution c_idx, + const std::vector& popold) const +{ + const auto &prob = pop.get_problem(); + auto dim = prob.get_nx(); + auto NP = pop.size(); + + std::vector r(5); // indexes of 5 selected population members + vector_double tmp(dim); + + /*-----We select at random 5 indexes from the population---------------------------------*/ + std::vector idxs(NP); + std::iota(idxs.begin(), idxs.end(), vector_double::size_type(0u)); + for (auto j = 0u; j < 5u; ++j) { // Durstenfeld's algorithm to select 5 indexes at random + auto idx = std::uniform_int_distribution(0u, NP - 1u - j)(m_e); + r[j] = idxs[idx]; + std::swap(idxs[idx], idxs[NP - 1u - j]); + } + + /*-------DE/best/1/exp--------------------------------------------------------------------*/ + /*-------The oldest DE variant but still not bad. However, we have found several---------*/ + /*-------optimization problems where misconvergence occurs.-------------------------------*/ + if (m_variant == 1u) { + tmp = popold[i]; + auto n = c_idx(m_e); + auto L = 0u; + do { + tmp[n] = gbIter[n] + m_F * (popold[r[1]][n] - popold[r[2]][n]); + n = (n + 1u) % dim; + ++L; + } while ((drng(m_e) < m_CR) && (L < dim)); + } + + /*-------DE/rand/1/exp-------------------------------------------------------------------*/ + /*-------This is one of my favourite strategies. It works especially well when the-------*/ + /*-------"gbIter[]"-schemes experience misconvergence. Try e.g. m_F=0.7 and m_CR=0.5---------*/ + /*-------as a first guess.---------------------------------------------------------------*/ + else if (m_variant == 2u) { + tmp = popold[i]; + auto n = c_idx(m_e); + decltype(dim) L = 0u; + do { + tmp[n] = popold[r[0]][n] + m_F * (popold[r[1]][n] - popold[r[2]][n]); + n = (n + 1u) % dim; + ++L; + } while ((drng(m_e) < m_CR) && (L < dim)); + } + /*-------DE/rand-to-best/1/exp-----------------------------------------------------------*/ + /*-------This variant seems to be one of the best strategies. Try m_F=0.85 and m_CR=1.------*/ + /*-------If you get misconvergence try to increase NP. If this doesn't help you----------*/ + /*-------should play around with all three control variables.----------------------------*/ + else if (m_variant == 3u) { + tmp = popold[i]; + auto n = c_idx(m_e); + auto L = 0u; + do { + tmp[n] = tmp[n] + m_F * (gbIter[n] - tmp[n]) + m_F * (popold[r[0]][n] - popold[r[1]][n]); + n = (n + 1u) % dim; + ++L; + } while ((drng(m_e) < m_CR) && (L < dim)); + } + /*-------DE/best/2/exp is another powerful variant worth trying--------------------------*/ + else if (m_variant == 4u) { + tmp = popold[i]; + auto n = c_idx(m_e); + auto L = 0u; + do { + tmp[n] = gbIter[n] + (popold[r[0]][n] + popold[r[1]][n] - popold[r[2]][n] - popold[r[3]][n]) * m_F; + n = (n + 1u) % dim; + ++L; + } while ((drng(m_e) < m_CR) && (L < dim)); + } + /*-------DE/rand/2/exp seems to be a robust optimizer for many functions-------------------*/ + else if (m_variant == 5u) { + tmp = popold[i]; + auto n = c_idx(m_e); + auto L = 0u; + do { + tmp[n] = popold[r[4]][n] + (popold[r[0]][n] + popold[r[1]][n] - popold[r[2]][n] - popold[r[3]][n]) * m_F; + n = (n + 1u) % dim; + ++L; + } while ((drng(m_e) < m_CR) && (L < dim)); + } + + /*=======Essentially same strategies but BINOMIAL CROSSOVER===============================*/ + /*-------DE/best/1/bin--------------------------------------------------------------------*/ + else if (m_variant == 6u) { + tmp = popold[i]; + auto n = c_idx(m_e); + for (decltype(dim) L = 0u; L < dim; ++L) { /* perform Dc binomial trials */ + if ((drng(m_e) < m_CR) || L + 1u == dim) { /* change at least one parameter */ + tmp[n] = gbIter[n] + m_F * (popold[r[1]][n] - popold[r[2]][n]); + } + n = (n + 1u) % dim; + } + } + /*-------DE/rand/1/bin-------------------------------------------------------------------*/ + else if (m_variant == 7u) { + tmp = popold[i]; + auto n = c_idx(m_e); + for (decltype(dim) L = 0u; L < dim; ++L) { /* perform Dc binomial trials */ + if ((drng(m_e) < m_CR) || L + 1u == dim) { /* change at least one parameter */ + tmp[n] = popold[r[0]][n] + m_F * (popold[r[1]][n] - popold[r[2]][n]); + } + n = (n + 1u) % dim; + } + } + /*-------DE/rand-to-best/1/bin-----------------------------------------------------------*/ + else if (m_variant == 8u) { + tmp = popold[i]; + auto n = c_idx(m_e); + for (decltype(dim) L = 0u; L < dim; ++L) { /* perform Dc binomial trials */ + if ((drng(m_e) < m_CR) || L + 1u == dim) { /* change at least one parameter */ + tmp[n] = tmp[n] + m_F * (gbIter[n] - tmp[n]) + m_F * (popold[r[0]][n] - popold[r[1]][n]); + } + n = (n + 1u) % dim; + } + } + /*-------DE/best/2/bin--------------------------------------------------------------------*/ + else if (m_variant == 9u) { + tmp = popold[i]; + auto n = c_idx(m_e); + for (decltype(dim) L = 0u; L < dim; ++L) { /* perform Dc binomial trials */ + if ((drng(m_e) < m_CR) || L + 1u == dim) { /* change at least one parameter */ + tmp[n] = gbIter[n] + (popold[r[0]][n] + popold[r[1]][n] - popold[r[2]][n] - popold[r[3]][n]) * m_F; + } + n = (n + 1u) % dim; + } + } + /*-------DE/rand/2/bin--------------------------------------------------------------------*/ + else if (m_variant == 10u) { + tmp = popold[i]; + auto n = c_idx(m_e); + for (decltype(dim) L = 0u; L < dim; ++L) { /* perform Dc binomial trials */ + if ((drng(m_e) < m_CR) || L + 1u == dim) { /* change at least one parameter */ + tmp[n] + = popold[r[4]][n] + (popold[r[0]][n] + popold[r[1]][n] - popold[r[2]][n] - popold[r[3]][n]) * m_F; + } + n = (n + 1u) % dim; + } + } + + return tmp; +} + +void de::update_pop(population &pop, const vector_double &newfitness, population::size_type i, + std::vector &fit, vector_double &gbfit, vector_double &gbX, + std::vector &popnew, const std::vector &popold, + const vector_double &tmp) const +{ + if (newfitness[0] <= fit[i][0]) { /* improved objective function value ? */ + fit[i] = newfitness; + popnew[i] = tmp; + // updates the individual in pop (avoiding to recompute the objective function) + pop.set_xf(i, popnew[i], newfitness); + + if (newfitness[0] <= gbfit[0]) { + /* if so...*/ + gbfit = newfitness; /* reset gbfit to new low...*/ + gbX = popnew[i]; + } + } else { + popnew[i] = popold[i]; + } +} + /// Sets the seed /** * @param seed the seed controlling the algorithm stochastic behaviour @@ -458,7 +459,7 @@ std::string de::get_extra_info() const template void de::serialize(Archive &ar, unsigned) { - detail::archive(ar, m_gen, m_F, m_CR, m_variant, m_Ftol, m_xtol, m_e, m_seed, m_verbosity, m_log); + detail::archive(ar, m_gen, m_F, m_CR, m_variant, m_Ftol, m_xtol, m_e, m_seed, m_verbosity, m_log, m_bfe); } } // namespace pagmo diff --git a/tests/de.cpp b/tests/de.cpp index 8ddeb566c..b625d7bbf 100644 --- a/tests/de.cpp +++ b/tests/de.cpp @@ -176,34 +176,18 @@ BOOST_AUTO_TEST_CASE(de_serialization_test) } } - BOOST_AUTO_TEST_CASE(bfe_usage_test) { - - std::chrono::steady_clock::time_point begin_1 = std::chrono::steady_clock::now(); - population pop{rosenbrock{10u}, 40u, 23u}; + population pop{rosenbrock{10u}, 200u, 23u}; de user_algo{1000000u, 0.7, 0.5, 2u, 1e-3, 1e-50, 23u}; user_algo.set_verbosity(1u); - user_algo.set_seed(23u); - user_algo.set_bfe(bfe{}); + user_algo.set_bfe(bfe{}); // This will use the default bfe. pop = user_algo.evolve(pop); - std::chrono::steady_clock::time_point end_1 = std::chrono::steady_clock::now(); - BOOST_CHECK(user_algo.get_log().size() < 6000u); - std::chrono::steady_clock::time_point begin_2 = std::chrono::steady_clock::now(); - population pop_2{rosenbrock{10u}, 40u, 23u}; + population pop_2{rosenbrock{10u}, 200u, 23u}; de user_algo_2{1000000u, 0.7, 0.5, 2u, 1e-3, 1e-50, 23u}; user_algo_2.set_verbosity(1u); pop_2 = user_algo_2.evolve(pop_2); - std::chrono::steady_clock::time_point end_2 = std::chrono::steady_clock::now(); - BOOST_CHECK(user_algo_2.get_log().size() < 6000u); - - BOOST_CHECK_EQUAL(pop.champion_f()[0], pop_2.champion_f()[0]); - - - std::cout << "Time difference = " << std::chrono::duration_cast(end_1 - begin_1).count() - << "[ms]" << std::endl; - std::cout << "Time difference = " << std::chrono::duration_cast(end_2 - begin_2).count() - << "[ms]" << std::endl; + BOOST_CHECK_EQUAL(pop.champion_f()[0], pop_2.champion_f()[0]); } From 6c91ffcfd6c01066306bd7e05cceca5d7f97096b Mon Sep 17 00:00:00 2001 From: Evgenii Vikulov Date: Tue, 6 Jun 2023 12:04:02 +0900 Subject: [PATCH 4/5] Code foramtting --- include/pagmo/algorithms/de.hpp | 4 ++-- src/algorithms/de.cpp | 12 ++++++------ tests/de.cpp | 3 +-- 3 files changed, 9 insertions(+), 10 deletions(-) diff --git a/include/pagmo/algorithms/de.hpp b/include/pagmo/algorithms/de.hpp index 86614b7e1..075ebd88e 100644 --- a/include/pagmo/algorithms/de.hpp +++ b/include/pagmo/algorithms/de.hpp @@ -163,7 +163,7 @@ class PAGMO_DLL_PUBLIC de return m_gen; } - // Sets the bfe + // Sets the bfe void set_bfe(const bfe &b); /// Algorithm name @@ -194,7 +194,7 @@ class PAGMO_DLL_PUBLIC de vector_double mutate(const population &pop, population::size_type i, const vector_double &gbIter, std::uniform_real_distribution drng, std::uniform_int_distribution c_idx, - const std::vector& popold) const; + const std::vector &popold) const; void update_pop(population &pop, const vector_double &newfitness, population::size_type i, std::vector &fit, vector_double &gbfit, vector_double &gbX, diff --git a/src/algorithms/de.cpp b/src/algorithms/de.cpp index 9ec7e6c1a..bf02c709b 100644 --- a/src/algorithms/de.cpp +++ b/src/algorithms/de.cpp @@ -137,7 +137,7 @@ population de::evolve(population pop) const // Main DE iterations for (decltype(m_gen) gen = 1u; gen <= m_gen; ++gen) { - + if (m_bfe) { // bfe is available: vector_double trial(NP * dim); @@ -158,7 +158,7 @@ population de::evolve(population pop) const ++pos; } } - + vector_double tmp(dim); vector_double newfitness(prob_f_dimension); auto fitnesses = (*m_bfe)(prob, trial); @@ -180,7 +180,7 @@ population de::evolve(population pop) const } } else { - + for (decltype(NP) i = 0u; i < NP; ++i) { auto tmp = mutate(pop, i, gbIter, drng, c_idx, popold); @@ -257,12 +257,12 @@ population de::evolve(population pop) const vector_double de::mutate(const population &pop, population::size_type i, const vector_double &gbIter, std::uniform_real_distribution drng, std::uniform_int_distribution c_idx, - const std::vector& popold) const + const std::vector &popold) const { - const auto &prob = pop.get_problem(); + const auto &prob = pop.get_problem(); auto dim = prob.get_nx(); auto NP = pop.size(); - + std::vector r(5); // indexes of 5 selected population members vector_double tmp(dim); diff --git a/tests/de.cpp b/tests/de.cpp index b625d7bbf..042aea4a6 100644 --- a/tests/de.cpp +++ b/tests/de.cpp @@ -32,7 +32,6 @@ see https://www.gnu.org/licenses/. */ #include #include -#include #include #include @@ -188,6 +187,6 @@ BOOST_AUTO_TEST_CASE(bfe_usage_test) de user_algo_2{1000000u, 0.7, 0.5, 2u, 1e-3, 1e-50, 23u}; user_algo_2.set_verbosity(1u); pop_2 = user_algo_2.evolve(pop_2); - + BOOST_CHECK_EQUAL(pop.champion_f()[0], pop_2.champion_f()[0]); } From 1761efc4f8d8a0c6e5f43bc9effbd7cc3a7234b0 Mon Sep 17 00:00:00 2001 From: Evgenii Vikulov Date: Tue, 6 Jun 2023 16:53:00 +0900 Subject: [PATCH 5/5] Turn off tests --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index e546841bf..2364363b4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -22,7 +22,7 @@ message(STATUS "System name: ${CMAKE_SYSTEM_NAME}") enable_testing() # Build option: enable the test suite. -option(PAGMO_BUILD_TESTS "Build the test suite." ON) +option(PAGMO_BUILD_TESTS "Build the test suite." OFF) # Build option: enable the benchmarks. option(PAGMO_BUILD_BENCHMARKS "Build the benchmarks." OFF)