Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fb fix pso memory #487

Open
wants to merge 8 commits into
base: master
Choose a base branch
from
Open

Fb fix pso memory #487

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
d9eb8ee
Keeping PSO variables at memory so that consecutive calls to evolve a…
Oct 5, 2021
e3b6d08
Refactoring: Move memory data to specific struct
Oct 5, 2021
003f435
Add PSO memory object serialization
Oct 6, 2021
e74abb3
Code Review: Replace population size type
Oct 6, 2021
c4ef9bf
Code Review: Header cleanup and build fixes
Oct 7, 2021
091c1f4
Code Review: Move m_V to memory
Oct 7, 2021
462e1c5
Replace broken paper link for genetic operators
Oct 7, 2021
9daa54f
Merge branch 'master' into fb-fix-pso-memory
darioizzo Apr 22, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 17 additions & 2 deletions include/pagmo/algorithms/pso.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ see https://www.gnu.org/licenses/. */
#ifndef PAGMO_ALGORITHMS_PSO_HPP
#define PAGMO_ALGORITHMS_PSO_HPP

#include <boost/optional.hpp>
#include <string>
#include <tuple>
#include <vector>
Expand Down Expand Up @@ -223,6 +224,21 @@ class PAGMO_DLL_PUBLIC pso
}

private:
struct memory {
std::vector<vector_double> m_V;
std::vector<vector_double> m_X;
std::vector<vector_double> m_lbX;
std::vector<vector_double> m_fit;
std::vector<vector_double> m_lbfit;
vector_double m_best_fit;
std::vector<std::vector<population::size_type>> m_neighb;
vector_double m_best_neighb;

// Object serialization
template <typename Archive>
void serialize(Archive &, unsigned);
};

// Object serialization
friend class boost::serialization::access;
template <typename Archive>
Expand Down Expand Up @@ -254,8 +270,7 @@ class PAGMO_DLL_PUBLIC pso
unsigned m_neighb_param;
// memory
bool m_memory;
// paricles' velocities
mutable std::vector<vector_double> m_V;
mutable boost::optional<pso::memory> m_memory_data;

mutable detail::random_engine_type m_e;
unsigned m_seed;
Expand Down
135 changes: 85 additions & 50 deletions src/algorithms/pso.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,10 @@ see https://www.gnu.org/licenses/. */
#include <pagmo/types.hpp>
#include <pagmo/utils/generic.hpp>

// NOTE: apparently this must be included *after*
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I find out that this is a "common problem" inside pagmo. This include block (comment+include) is copy-pasted in several algorithms. In order to maintain the pattern, I decided to copy-paste it here too.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The Boost.serialisation library is unfortunately sensitive to the header order inclusion when it comes to the serialisation of polymorphic classes via base pointers (which we need in the implementation of pagmo's type-erased wrappers).

// the other serialization headers.
#include <boost/serialization/optional.hpp>

namespace pagmo
{

Expand Down Expand Up @@ -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.) {
Expand Down Expand Up @@ -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<vector_double>(swarm_size, dummy);
std::vector<vector_double> 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;
tarcisiofischer marked this conversation as resolved.
Show resolved Hide resolved
double sum_forces = 0.;

double r1 = 0.;
double r2 = 0.;
Expand All @@ -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]);
}
}

Expand All @@ -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]);
}
}

Expand All @@ -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]);
}
}

Expand All @@ -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]);
}
}

Expand All @@ -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]));
}
}

Expand All @@ -319,36 +341,36 @@ 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<double>(neighb[p].size()));
V[p][d] = m_omega * (V[p][d] + sum_forces / static_cast<double>(neighb[p].size()));
}
}

// We now check that the velocity does not exceed the maximum allowed per component
// 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<double>(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<double>(lb[d],ub[d])(m_drng);
// V[p][d] = new_x - X[p][d];
}
Expand Down Expand Up @@ -393,13 +415,13 @@ population pso::evolve(population pop) const
auto best = local_fits[static_cast<vector_double::size_type>(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<double>(m_V[i].size());
mean_velocity /= static_cast<double>(V[i].size());
}
// We compute the average distance across particles (NOTE: N^2 complexity)
auto avg_dist = 0.;
Expand Down Expand Up @@ -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;
}

Expand Down Expand Up @@ -481,8 +509,15 @@ std::string pso::get_extra_info() const
template <typename Archive>
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 <typename Archive>
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);
}

/**
Expand Down
Loading