diff --git a/include/pagmo/algorithms/pso.hpp b/include/pagmo/algorithms/pso.hpp index c36a15b75..11861c0f8 100644 --- a/include/pagmo/algorithms/pso.hpp +++ b/include/pagmo/algorithms/pso.hpp @@ -29,6 +29,7 @@ see https://www.gnu.org/licenses/. */ #ifndef PAGMO_ALGORITHMS_PSO_HPP #define PAGMO_ALGORITHMS_PSO_HPP +#include #include #include #include @@ -223,6 +224,21 @@ class PAGMO_DLL_PUBLIC pso } private: + struct memory { + std::vector m_V; + std::vector m_X; + std::vector m_lbX; + std::vector m_fit; + std::vector m_lbfit; + vector_double m_best_fit; + std::vector> m_neighb; + vector_double m_best_neighb; + + // Object serialization + template + void serialize(Archive &, unsigned); + }; + // Object serialization friend class boost::serialization::access; template @@ -254,8 +270,7 @@ class PAGMO_DLL_PUBLIC pso unsigned m_neighb_param; // memory bool m_memory; - // paricles' velocities - mutable std::vector m_V; + mutable boost::optional m_memory_data; mutable detail::random_engine_type m_e; unsigned m_seed; diff --git a/src/algorithms/pso.cpp b/src/algorithms/pso.cpp index a8ba39d7a..7cc487836 100644 --- a/src/algorithms/pso.cpp +++ b/src/algorithms/pso.cpp @@ -46,6 +46,10 @@ see https://www.gnu.org/licenses/. */ #include #include +// NOTE: apparently this must be included *after* +// the other serialization headers. +#include + namespace pagmo { @@ -73,7 +77,7 @@ constexpr int pso_vonNeumann_neighb_diff[4][2] = {{-1, 0}, {1, 0}, {0, -1}, {0, pso::pso(unsigned gen, double omega, double eta1, double eta2, double max_vel, unsigned variant, unsigned neighb_type, unsigned neighb_param, bool memory, unsigned seed) : m_max_gen(gen), m_omega(omega), m_eta1(eta1), m_eta2(eta2), m_max_vel(max_vel), m_variant(variant), - m_neighb_type(neighb_type), m_neighb_param(neighb_param), m_memory(memory), m_V(), m_e(seed), m_seed(seed), + m_neighb_type(neighb_type), m_neighb_param(neighb_param), m_memory(memory), m_e(seed), m_seed(seed), m_verbosity(0u), m_log() { if (m_omega < 0. || m_omega > 1.) { @@ -179,44 +183,62 @@ population pso::evolve(population pop) const } // Copy the particle positions and their fitness - for (decltype(swarm_size) i = 0u; i < swarm_size; ++i) { - X[i] = pop.get_x()[i]; - lbX[i] = pop.get_x()[i]; + // If calling from memory, the positions from last run may not be the same as the best population + // so it is make a correction here + if (m_memory && m_memory_data) { + X = m_memory_data->m_X; + lbX = m_memory_data->m_lbX; + + fit = m_memory_data->m_fit; + lbfit = m_memory_data->m_lbfit; + } else { + for (decltype(swarm_size) i = 0u; i < swarm_size; ++i) { + X[i] = pop.get_x()[i]; + lbX[i] = pop.get_x()[i]; - fit[i] = pop.get_f()[i]; - lbfit[i] = pop.get_f()[i]; + fit[i] = pop.get_f()[i]; + lbfit[i] = pop.get_f()[i]; + } } // Initialize the particle velocities if necessary - if ((m_V.size() != swarm_size) || (!m_memory)) { - m_V = std::vector(swarm_size, dummy); + std::vector V(swarm_size, dummy); + if (m_memory && m_memory_data) { + V = m_memory_data->m_V; + } else { for (decltype(swarm_size) i = 0u; i < swarm_size; ++i) { for (decltype(dim) j = 0u; j < dim; ++j) { - m_V[i][j] = uniform_real_from_range(minv[j], maxv[j], m_e); + V[i][j] = uniform_real_from_range(minv[j], maxv[j], m_e); } } } // Initialize the Swarm's topology - switch (m_neighb_type) { - case 1: - initialize_topology__gbest(pop, best_neighb, best_fit, neighb); - break; - case 3: - initialize_topology__von(neighb); - break; - case 4: - initialize_topology__adaptive_random(neighb); - // need to track improvements in best found fitness, to know when to rewire - best_fit = pop.get_f()[pop.best_idx()]; - break; - case 2: - default: - initialize_topology__lbest(neighb); + if (m_memory && m_memory_data) { + neighb = m_memory_data->m_neighb; + best_fit = m_memory_data->m_best_fit; + best_neighb = m_memory_data->m_best_neighb; + } else { + switch (m_neighb_type) { + case 1: + initialize_topology__gbest(pop, best_neighb, best_fit, neighb); + break; + case 3: + initialize_topology__von(neighb); + break; + case 4: + initialize_topology__adaptive_random(neighb); + // need to track improvements in best found fitness, to know when to rewire + best_fit = pop.get_f()[pop.best_idx()]; + break; + case 2: + default: + initialize_topology__lbest(neighb); + } } // auxiliary variables specific to the Fully Informed Particle Swarm variant double acceleration_coefficient = m_eta1 + m_eta2; - double sum_forces; + double sum_forces = 0.; double r1 = 0.; double r2 = 0.; @@ -242,8 +264,8 @@ population pso::evolve(population pop) const for (decltype(dim) d = 0u; d < dim; ++d) { r1 = drng(m_e); r2 = drng(m_e); - m_V[p][d] = m_omega * m_V[p][d] + m_eta1 * r1 * (lbX[p][d] - X[p][d]) - + m_eta2 * r2 * (best_neighb[d] - X[p][d]); + V[p][d] = m_omega * V[p][d] + m_eta1 * r1 * (lbX[p][d] - X[p][d]) + + m_eta2 * r2 * (best_neighb[d] - X[p][d]); } } @@ -253,8 +275,8 @@ population pso::evolve(population pop) const else if (m_variant == 2u) { for (decltype(dim) d = 0u; d < dim; ++d) { r1 = drng(m_e); - m_V[p][d] = m_omega * m_V[p][d] + m_eta1 * r1 * (lbX[p][d] - X[p][d]) - + m_eta2 * r1 * (best_neighb[d] - X[p][d]); + V[p][d] = m_omega * V[p][d] + m_eta1 * r1 * (lbX[p][d] - X[p][d]) + + m_eta2 * r1 * (best_neighb[d] - X[p][d]); } } @@ -264,8 +286,8 @@ population pso::evolve(population pop) const r1 = drng(m_e); r2 = drng(m_e); for (decltype(dim) d = 0u; d < dim; ++d) { - m_V[p][d] = m_omega * m_V[p][d] + m_eta1 * r1 * (lbX[p][d] - X[p][d]) - + m_eta2 * r2 * (best_neighb[d] - X[p][d]); + V[p][d] = m_omega * V[p][d] + m_eta1 * r1 * (lbX[p][d] - X[p][d]) + + m_eta2 * r2 * (best_neighb[d] - X[p][d]); } } @@ -275,8 +297,8 @@ population pso::evolve(population pop) const else if (m_variant == 4u) { r1 = drng(m_e); for (decltype(dim) d = 0u; d < dim; ++d) { - m_V[p][d] = m_omega * m_V[p][d] + m_eta1 * r1 * (lbX[p][d] - X[p][d]) - + m_eta2 * r1 * (best_neighb[d] - X[p][d]); + V[p][d] = m_omega * V[p][d] + m_eta1 * r1 * (lbX[p][d] - X[p][d]) + + m_eta2 * r1 * (best_neighb[d] - X[p][d]); } } @@ -297,9 +319,9 @@ population pso::evolve(population pop) const for (decltype(dim) d = 0u; d < dim; ++d) { r1 = drng(m_e); r2 = drng(m_e); - m_V[p][d] = m_omega - * (m_V[p][d] + m_eta1 * r1 * (lbX[p][d] - X[p][d]) - + m_eta2 * r2 * (best_neighb[d] - X[p][d])); + V[p][d] + = m_omega + * (V[p][d] + m_eta1 * r1 * (lbX[p][d] - X[p][d]) + m_eta2 * r2 * (best_neighb[d] - X[p][d])); } } @@ -319,7 +341,7 @@ population pso::evolve(population pop) const for (decltype(neighb[p].size()) n = 0u; n < neighb[p].size(); ++n) { sum_forces += drng(m_e) * acceleration_coefficient * (lbX[neighb[p][n]][d] - X[p][d]); } - m_V[p][d] = m_omega * (m_V[p][d] + sum_forces / static_cast(neighb[p].size())); + V[p][d] = m_omega * (V[p][d] + sum_forces / static_cast(neighb[p].size())); } } @@ -327,28 +349,28 @@ population pso::evolve(population pop) const // and we perform the position update and the feasibility correction for (decltype(dim) d = 0u; d < dim; ++d) { - if (m_V[p][d] > maxv[d]) { - m_V[p][d] = maxv[d]; + if (V[p][d] > maxv[d]) { + V[p][d] = maxv[d]; } - else if (m_V[p][d] < minv[d]) { - m_V[p][d] = minv[d]; + else if (V[p][d] < minv[d]) { + V[p][d] = minv[d]; } // update position - new_x = X[p][d] + m_V[p][d]; + new_x = X[p][d] + V[p][d]; // feasibility correction // (velocity updated to that which would have taken the previous position // to the newly corrected feasible position) if (new_x < lb[d]) { new_x = lb[d]; - m_V[p][d] = 0.; + V[p][d] = 0.; // new_x = boost::uniform_real(lb[d],ub[d])(m_drng); // V[p][d] = new_x - X[p][d]; } else if (new_x > ub[d]) { new_x = ub[d]; - m_V[p][d] = 0.; + V[p][d] = 0.; // new_x = boost::uniform_real(lb[d],ub[d])(m_drng); // V[p][d] = new_x - X[p][d]; } @@ -393,13 +415,13 @@ population pso::evolve(population pop) const auto best = local_fits[static_cast(idx_best)]; // We compute a measure for the average particle velocity across the swarm auto mean_velocity = 0.; - for (decltype(m_V.size()) i = 0u; i < m_V.size(); ++i) { - for (decltype(m_V[i].size()) j = 0u; j < m_V[i].size(); ++j) { + for (decltype(V.size()) i = 0u; i < V.size(); ++i) { + for (decltype(V[i].size()) j = 0u; j < V[i].size(); ++j) { if (ub[j] > lb[j]) { - mean_velocity += std::abs(m_V[i][j] / (ub[j] - lb[j])); + mean_velocity += std::abs(V[i][j] / (ub[j] - lb[j])); } // else 0 } - mean_velocity /= static_cast(m_V[i].size()); + mean_velocity /= static_cast(V[i].size()); } // We compute the average distance across particles (NOTE: N^2 complexity) auto avg_dist = 0.; @@ -439,6 +461,12 @@ population pso::evolve(population pop) const for (decltype(swarm_size) i = 0u; i < swarm_size; ++i) { pop.set_xf(i, lbX[i], lbfit[i]); } + + // Keep memory variables only if asked for + if (m_memory) { + m_memory_data = pso::memory{V, X, lbX, fit, lbfit, best_fit, neighb, best_neighb}; + } + return pop; } @@ -481,8 +509,15 @@ std::string pso::get_extra_info() const template void pso::serialize(Archive &ar, unsigned) { - detail::archive(ar, m_max_gen, m_omega, m_eta1, m_eta2, m_max_vel, m_variant, m_neighb_type, m_neighb_param, m_e, - m_seed, m_verbosity, m_log); + detail::archive(ar, m_max_gen, m_omega, m_eta1, m_eta2, m_max_vel, m_variant, m_neighb_type, m_neighb_param, + m_memory, m_memory_data, m_e, m_seed, m_verbosity, m_log); +} + +// Object'm memory serialization +template +void pso::memory::serialize(Archive &ar, unsigned) +{ + detail::archive(ar, m_V, m_X, m_lbX, m_fit, m_lbfit, m_best_fit, m_neighb, m_best_neighb); } /** diff --git a/tests/pso.cpp b/tests/pso.cpp index 122111047..384da8c03 100644 --- a/tests/pso.cpp +++ b/tests/pso.cpp @@ -134,6 +134,7 @@ BOOST_AUTO_TEST_CASE(evolve_test) } } } + BOOST_AUTO_TEST_CASE(setters_getters_test) { pso user_algo{5000u, 0.79, 2., 2., 0.1, 5u, 2u, 4u, false, 23u}; @@ -184,3 +185,99 @@ BOOST_AUTO_TEST_CASE(serialization_test) BOOST_CHECK_CLOSE(std::get<5>(before_log[i]), std::get<5>(after_log[i]), 1e-8); } } + +BOOST_AUTO_TEST_CASE(pso_memory_test) +{ + // We check here that when memory is true calling evolve(pop) two times on 1 gen + // is the same as calling 1 time evolve with 2 gens + auto omega = 0.5; + auto eta1 = 0.5; + auto eta2 = 0.5; + auto max_vel = 0.5; + auto neighb_param = 4u; + auto seed = 42u; + auto pop_size = 10u; + + for (auto variant = 1u; variant <= 6; ++variant) { + for (auto neighb_type = 1u; neighb_type <= 4; ++neighb_type) { + auto log1 = ([&]() { + auto n_generations = 1u; + auto memory = true; + algorithm algo{ + pso{n_generations, omega, eta1, eta2, max_vel, variant, neighb_type, neighb_param, memory, seed}}; + algo.set_verbosity(1u); + problem prob{rosenbrock{25u}}; + population pop{prob, pop_size, seed}; + pop = algo.evolve(pop); + pop = algo.evolve(pop); + pop = algo.evolve(pop); + return algo.extract()->get_log(); + })(); + + auto log2 = ([&]() { + auto n_generations = 3u; + auto memory = false; + algorithm algo{ + pso{n_generations, omega, eta1, eta2, max_vel, variant, neighb_type, neighb_param, memory, seed}}; + algo.set_verbosity(1u); + problem prob{rosenbrock{25u}}; + population pop{prob, pop_size, seed}; + pop = algo.evolve(pop); + return algo.extract()->get_log(); + })(); + + // indexes 0 (Gen) and 1 (Fevals) may be different + // Check index 2 (gbest) + BOOST_CHECK_CLOSE(std::get<2>(log1[0]), std::get<2>(log2[2]), 1e-8); + // Check index 3 (Mean Vel.) + BOOST_CHECK_CLOSE(std::get<3>(log1[0]), std::get<3>(log2[2]), 1e-8); + // Check index 4 (Mean lbest) + BOOST_CHECK_CLOSE(std::get<4>(log1[0]), std::get<4>(log2[2]), 1e-8); + // Check index 5 (Avg. Dist.) + BOOST_CHECK_CLOSE(std::get<5>(log1[0]), std::get<5>(log2[2]), 1e-8); + + // Make sure serializing and deserializing the results will give the same answer + auto serialize_and_deserialize = [](algorithm &algo) { + std::stringstream ss; + { + boost::archive::binary_oarchive oarchive(ss); + oarchive << algo; + } + algo = algorithm{}; + { + boost::archive::binary_iarchive iarchive(ss); + iarchive >> algo; + } + }; + auto log3 = ([&]() { + auto n_generations = 1u; + auto memory = true; + algorithm algo{ + pso{n_generations, omega, eta1, eta2, max_vel, variant, neighb_type, neighb_param, memory, seed}}; + algo.set_verbosity(1u); + problem prob{rosenbrock{25u}}; + population pop{prob, pop_size, seed}; + pop = algo.evolve(pop); + serialize_and_deserialize(algo); + pop = algo.evolve(pop); + serialize_and_deserialize(algo); + pop = algo.evolve(pop); + return algo.extract()->get_log(); + })(); + // While comparing the serialized vs non-serialized versions, all logs must be the same. + BOOST_CHECK_EQUAL(log1.size(), log3.size()); + // Check index 0 (Gen) + BOOST_CHECK_EQUAL(std::get<0>(log1[0]), std::get<0>(log3[0])); + // Check index 1 (Fevals) + BOOST_CHECK_EQUAL(std::get<1>(log1[0]), std::get<1>(log3[0])); + // Check index 2 (gbest) + BOOST_CHECK_CLOSE(std::get<2>(log1[0]), std::get<2>(log3[0]), 1e-8); + // Check index 3 (Mean Vel.) + BOOST_CHECK_CLOSE(std::get<3>(log1[0]), std::get<3>(log3[0]), 1e-8); + // Check index 4 (Mean lbest) + BOOST_CHECK_CLOSE(std::get<4>(log1[0]), std::get<4>(log3[0]), 1e-8); + // Check index 5 (Avg. Dist.) + BOOST_CHECK_CLOSE(std::get<5>(log1[0]), std::get<5>(log3[0]), 1e-8); + } + } +}