Fb fix pso memory

include/pagmo/algorithms/pso.hpp

@@ -31,6 +31,7 @@ see https://www.gnu.org/licenses/. */
 
 #include <string>
 #include <tuple>
+#include <type_traits>
 #include <vector>
 
 #include <pagmo/algorithm.hpp>
@@ -223,6 +224,18 @@ class PAGMO_DLL_PUBLIC pso
     }
 
 private:
+    struct memory {
+        using population_size_t = decltype(std::declval<population>().size());
+
+        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_t>> m_neighb;
+        vector_double m_best_neighb;
+    };
+
     // Object serialization
     friend class boost::serialization::access;
     template <typename Archive>
@@ -254,9 +267,10 @@ class PAGMO_DLL_PUBLIC pso
     unsigned m_neighb_param;
     // memory
     bool m_memory;
+    mutable std::optional<pso::memory> m_memory_data;
+
     // paricles' velocities
     mutable std::vector<vector_double> m_V;
-
     mutable detail::random_engine_type m_e;
     unsigned m_seed;
     unsigned m_verbosity;

src/algorithms/pso.cpp

@@ -179,12 +179,22 @@ 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
@@ -198,25 +208,31 @@ population pso::evolve(population pop) const
     }
 
     // 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.;
@@ -439,6 +455,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{X, lbX, fit, lbfit, best_fit, neighb, best_neighb};
+    }
+
     return pop;
 }
 

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,54 @@ 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.5d;
+    auto eta1 = 0.5d;
+    auto eta2 = 0.5d;
+    auto max_vel = 0.5d;
+    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;
+                pso user_algo{n_generations, omega, eta1, eta2, max_vel, variant, neighb_type, neighb_param, memory, seed};
+                user_algo.set_verbosity(1u);
+                problem prob{rosenbrock{25u}};
+                population pop{prob, pop_size, seed};
+                pop = user_algo.evolve(pop);
+                pop = user_algo.evolve(pop);
+                pop = user_algo.evolve(pop);
+                return user_algo.get_log();
+            })();
+
+            auto log2 = ([&]() {
+                auto n_generations = 3u;
+                auto memory = false;
+                pso user_algo{n_generations, omega, eta1, eta2, max_vel, variant, neighb_type, neighb_param, memory, seed};
+                user_algo.set_verbosity(1u);
+                problem prob{rosenbrock{25u}};
+                population pop{prob, pop_size, seed};
+                pop = user_algo.evolve(pop);
+                return user_algo.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);
+        }
+    }
+}