Rough working implementation of parallel particle swarm optimization using MPI

CMakeLists.txt

@@ -2,14 +2,19 @@ cmake_minimum_required(VERSION 3.2)
 
 project(Oersted)
 
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=gnu++17")
-#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O3 --coverage -fopenmp")
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O3 -fopenmp")
+find_package(MPI REQUIRED)
+find_package(OpenMP)
 
+include_directories(${MPI_INCLUDE_PATH})
 include_directories(./lib/)
 include_directories(./lib/Eigen/)
 include_directories(./lib/GoogleTest/googletest/include/gtest/)
 
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=gnu++17")
+#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O3 --coverage -fopenmp")
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O3 -fopenmp")
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+
 add_subdirectory(./lib/)
 
 include_directories(./src/)

src/Optimization/src/Particle.cpp

 #include "Particle.h"
 
-void Particle::update_personal_best(std::map<std::string,double_t> go) {
-        double_t value = LocalObjective(go);
-
-        if (value < PersonalBest) {
-            HasImproved = true;
-            PersonalBest = value;
-            for (auto &key_state : State) {
-                auto const& key = key_state.first;
-                personal_best(key) = position(key);
-            }
+void Particle::update_personal_best(ObjectiveMap go) {
+    update_personal_best(LocalObjective(go));
+}
+
+void Particle::update_personal_best(double_t value) {
+    if (value < PersonalBest) {
+        HasImproved = true;
+        PersonalBest = value;
+        for (auto &ki : State.map()) {
+            auto const& k = ki.first;
+            personal_best(k) = position(k);
         }
+    }
 }
 
-void Particle::update_global_best(std::map<std::string,double_t> go, Particle &p) {
-    double_t value = LocalObjective(go);
+void Particle::update_global_best(ObjectiveMap go, Particle &p) {
+    update_global_best(LocalObjective(go),p);
+}
 
+void Particle::update_global_best(double_t value, Particle &p) {
     if (value < GlobalBest) {
         HasImproved = true;
         GlobalBest = value;
-        for (auto &key_state : State) {
-            std::string const &key = key_state.first;
-            global_best(key) = p.position(key);
+        for (auto &ki : State.map()) {
+            auto const &k = ki.first;
+            global_best(k) = p.position(k);
         }
     }
 }
@@ -29,16 +33,14 @@ void Particle::update_global_best(std::map<std::string,double_t> go, Particle &p
 void Particle::update_velocity(std::mt19937 rng, double_t omega, double_t phip, double_t phig) {
     std::uniform_real_distribution<> dist(0.0,1.0);
 
-    for (auto &key_state : State) {
-        CoordinateState &state = key_state.second;
-        state.update_velocity(omega, dist(rng) * phip, dist(rng) * phig);
+    for (CoordinateState &s : State.vector()) {
+        s.update_velocity(omega, dist(rng) * phip, dist(rng) * phig);
     }
 }
 
 void Particle::update_position() {
-    for (auto &key_state : State) {
-        CoordinateState &state = key_state.second;
-        state.update_position();
+    for (CoordinateState &s : State.vector() ) {
+        s.update_position();
     }
 }
 
@@ -53,11 +55,11 @@ bool Particle::is_converged(CoordinateSpace &space, double_t state_tol, double_t
         return false;
     }
 
-    for (auto &key_state : State) {
-        std::string key = key_state.first;
-        CoordinateState &state = key_state.second;
+    for (auto &ki : State.map()) {
+        std::string k = ki.first;
+        CoordinateState &s = state(k);
 
-        if(!state.is_converged(state_tol * space.coordinate(key).range())) {
+        if(!s.is_converged(state_tol * space.coordinate(k).range())) {
             return false;
         }
     }
@@ -65,8 +67,12 @@ bool Particle::is_converged(CoordinateSpace &space, double_t state_tol, double_t
     return true;
 }
 
-std::ostream& operator<<(std::ostream& os, const Particle& p) {
-    os << "P-Best: " << p.PersonalBest << ", G-Best: " << p.GlobalBest;
+std::ostream& operator<<(std::ostream& os, Particle& p) {
+    for (auto &m : p.state().map()) {
+        os << m.first << ": " << p.state(m.second).PersonalBest << ", ";
+    }
+
+    os << std:: endl << "P-Best: " << p.PersonalBest << ", G-Best: " << p.GlobalBest;
 }
 
 Particle CoordinateSpace::new_particle(std::mt19937 &rng) {
@@ -137,7 +143,6 @@ void Swarm::run() {
 
     bool all_converged{false};
     for (size_t iter = 0; iter != MaximumIterations && !all_converged; ++iter) {
-        std::cout << iter << std::endl;
         for (size_t i = 0; i != Particles.size(); ++i) {
             auto go = GlobalObjective(Particles[i]);
 
@@ -172,4 +177,143 @@ void Swarm::run() {
             }
         }
     }
+}
+
+void Swarm::mpirun() {
+    int TAG_PARTICLE_INDEX{1};
+    int TAG_PARTICLE_DATA{2};
+    int TAG_OBJECTIVE_DATA{3};
+    int TAG_EXIT{4};
+    int MPI_RANK;
+    int MPI_SIZE;
+    int FLAG{0};
+
+    MPI_Datatype MPI_CS;
+    MPI_Type_contiguous(sizeof(CoordinateState) / sizeof(MPI_DOUBLE), MPI_DOUBLE, &MPI_CS);
+    MPI_Type_commit(&MPI_CS);
+
+    MPI_Datatype MPI_STATE;
+    MPI_Type_vector((int)Particles[0].state().size(),1,1,MPI_CS,&MPI_STATE);
+    MPI_Type_commit(&MPI_STATE);
+
+
+    MPI_Comm_size(MPI_COMM_WORLD, &MPI_SIZE);
+    MPI_Comm_rank(MPI_COMM_WORLD, &MPI_RANK);
+
+    MPI_Status status;
+    MPI_Request request;
+
+    std::vector<int> PARTICLE_NUMBERS(Particles.size());
+    std::iota(PARTICLE_NUMBERS.begin(),PARTICLE_NUMBERS.end(),0);
+
+    std::deque<int> queue(Particles.size());
+    std::iota(queue.begin(),queue.end(),0);
+
+    MPI_Barrier(MPI_COMM_WORLD);
+
+    if (MPI_RANK == 0) {
+        std::uniform_real_distribution<> dist(0.0,1.0);
+        std::vector<bool> mpi_ready(MPI_SIZE,true);
+        double_t iter{0.0};
+        bool all_ready{false};
+        bool all_converged{false};
+        while ((iter < MaximumIterations && !all_converged) || !all_ready) {
+            // Send particle data to mpi processes
+            for (int i = 1; i != MPI_SIZE && iter < MaximumIterations && !all_converged; ++i) {
+                if (mpi_ready[i] && !queue.empty()) {
+                    int p = queue.front();
+                    queue.pop_front();
+
+                    MPI_Send(PARTICLE_NUMBERS.data() + p, 1, MPI_INT, i, TAG_PARTICLE_INDEX, MPI_COMM_WORLD);
+
+                    MPI_Send(Particles[p].state().vector().data(), 1, MPI_STATE, i, TAG_PARTICLE_DATA, MPI_COMM_WORLD);
+
+                    mpi_ready[i] = false;
+                }
+            }
+
+            // Collect particle data from mpi processes
+            for (int i = 1; i != MPI_SIZE; ++i) {
+                MPI_Iprobe(i, TAG_PARTICLE_INDEX, MPI_COMM_WORLD, &FLAG, &status);
+                if (FLAG) {
+                    int q;
+                    std::vector<double_t> objective(Particles.size(), 0.0);
+
+                    MPI_Recv(&q, 1, MPI_INT, i, TAG_PARTICLE_INDEX, MPI_COMM_WORLD, &status);
+                    MPI_Recv(objective.data(), (int) objective.size(), MPI_DOUBLE, i, TAG_OBJECTIVE_DATA,
+                             MPI_COMM_WORLD, &status);
+
+                    // Update personal best
+                    for (size_t j = 0; j != Particles.size(); ++j) {
+                        Particles[j].update_global_best(objective[j], Particles[q]);
+                    }
+
+                    Particles[q].update_personal_best(objective[q]);
+                    if (!Particles[q].has_improved()) {
+                        Space.perturb(Particles[q],RNG);
+                    }
+
+                    // Update particle q velocity and position
+                    Omega = 2.0 / (1 + M_SQRT2) * dist(RNG);
+                    Particles[q].update_velocity(RNG, Omega, PhiP, PhiG);
+                    Space.limit_velocity(Particles[q]);
+                    Particles[q].update_position();
+                    Space.limit_position(Particles[q]);
+
+                    mpi_ready[i] = true;
+                    queue.push_back(q);
+
+                    iter = iter + 1.0 / Particles.size();
+                }
+            }
+
+            all_converged = true;
+            for (auto &p : Particles) {
+                all_converged &= p.is_converged(Space,CoordinateTolerance,ObjectiveTolerance);
+            }
+
+            all_ready = true;
+            for (bool b : mpi_ready) {
+                all_ready &= b;
+            }
+        }
+
+        // Exit loop
+        for (int i = 1; i != MPI_SIZE; ++i) {
+            MPI_Isend(NULL, 0, MPI_INT, i, TAG_EXIT, MPI_COMM_WORLD, &request);
+        }
+    }
+
+    while (MPI_RANK != 0) {
+        //int data;
+
+        MPI_Iprobe(0,TAG_EXIT,MPI_COMM_WORLD,&FLAG,&status);
+        if (FLAG) {
+            MPI_Irecv(NULL,0,MPI_INT,0,TAG_EXIT,MPI_COMM_WORLD,&request);
+            break;
+        }
+
+        MPI_Iprobe(0,TAG_PARTICLE_INDEX,MPI_COMM_WORLD,&FLAG,&status);
+        if (FLAG) {
+            int p;
+
+            MPI_Recv(&p,1,MPI_INT,0,TAG_PARTICLE_INDEX,MPI_COMM_WORLD,&status);
+            MPI_Recv(Particles[p].state().vector().data(),1,MPI_STATE,0,TAG_PARTICLE_DATA,MPI_COMM_WORLD,&status);
+
+            Particles[p].state(0).PersonalBest = MPI_RANK;
+
+            ObjectiveMap go = GlobalObjective(Particles[p]);
+
+            std::vector<double_t> objective(Particles.size(),0.0);
+            for (size_t i = 0; i != Particles.size(); ++i) {
+                objective[i] = Particles[i].local_objective(go);
+            }
+
+            MPI_Send(&p,1,MPI_INT,0,TAG_PARTICLE_INDEX,MPI_COMM_WORLD);
+            MPI_Send(objective.data(),(int)objective.size(),MPI_DOUBLE,0,TAG_OBJECTIVE_DATA,MPI_COMM_WORLD);
+        }
+    }
+
+    MPI_Type_free(&MPI_CS);
+    MPI_Type_free(&MPI_STATE);
 }
\ No newline at end of file

src/Optimization/src/Particle.h

@@ -8,12 +8,48 @@
 #include <functional>
 #include <random>
 #include <iostream>
+#include <deque>
+#include <mpi.h>
+
 
 class Particle;
 class CoordinateSpace;
 class Coordinate;
 class CoordinateState;
 
+template <class Key, class T>
+class ContiguousDataMap {
+public:
+    // TODO: Copy Constructor
+
+    T &operator[] (Key k) {
+        auto it = Map.find(k);
+        if (it != Map.end()) {
+            return Vector[it->second];
+        } else {
+            Map[k] = Vector.size();
+            Vector.push_back(T());
+            return Vector.back();
+        }
+    };
+
+    T &operator[] (size_t k) { return Vector[k]; };
+
+    std::map<Key,size_t> &map() { return Map; };
+
+    std::vector<T> &vector() { return Vector; };
+
+    size_t size() const { return Vector.size(); };
+
+protected:
+    std::map<Key,size_t> Map;
+
+    std::vector<T> Vector;
+};
+
+using ObjectiveMap = ContiguousDataMap<std::string,double_t>;
+using ParticleState = ContiguousDataMap<std::string,CoordinateState>;
+
 class CoordinateState {
 public:
     CoordinateState() {};
@@ -57,27 +93,35 @@ protected:
 
 class Particle {
 public:
+    ParticleState &state() { return State; };
+
     CoordinateState &state(std::string key) { return State[key]; };
 
-    double_t &position(std::string key) { return State[key].Position; };
+    CoordinateState &state(size_t index) { return State[index]; };
 
-    double_t &velocity(std::string key) { return State[key].Velocity; };
+    double_t &position(std::string key) { return state(key).Position; };
 
-    double_t &personal_best(std::string key) { return State[key].PersonalBest; };
+    double_t &velocity(std::string key) { return state(key).Velocity; };
 
-    double_t &global_best(std::string key) { return State[key].GlobalBest; };
+    double_t &personal_best(std::string key) { return state(key).PersonalBest; };
+
+    double_t &global_best(std::string key) { return state(key).GlobalBest; };
 
     double_t &personal_best() { return PersonalBest; };
 
     double_t &global_best() { return GlobalBest; };
 
-    std::function<double_t(std::map<std::string,double_t>)> &local_objective() { return LocalObjective; };
+    std::function<double_t(ObjectiveMap)> &local_objective() { return LocalObjective; };
+
+    double_t local_objective(ObjectiveMap &go) { return LocalObjective(go); };
 
-    double_t local_objective(std::map<std::string,double_t> &go) { return LocalObjective(go); };
+    void update_personal_best(ObjectiveMap objectives);
 
-    void update_personal_best(std::map<std::string,double_t> objectives);
+    void update_personal_best(double_t value);
 
-    void update_global_best(std::map<std::string,double_t> objectives, Particle &p);
+    void update_global_best(ObjectiveMap objectives, Particle &p);
+
+    void update_global_best(double_t value, Particle &p);
 
     void update_velocity(std::mt19937 rng, double_t omega, double_t phip, double_t phig);
 
@@ -87,16 +131,16 @@ public:
 
     bool is_converged(CoordinateSpace &ss, double_t coordinate_tol, double_t objective_tol);
 
-    friend std::ostream& operator<<(std::ostream& os, const Particle& p);
+    friend std::ostream& operator<<(std::ostream& os, Particle& p);
 
-    std::function<double_t(std::map<std::string,double_t>)> LocalObjective;
+    std::function<double_t(ObjectiveMap)> LocalObjective;
 
     double_t PersonalBest{DBL_MAX};
 
     double_t GlobalBest{DBL_MAX};
 
 protected:
-    std::map<std::string,CoordinateState> State;
+    ParticleState State;
 
     bool HasImproved{true};
 
@@ -121,12 +165,12 @@ public:
 
     double_t SpeedLimit{0.2};
 
-    double_t Perturbation{0.2e-2};
+    double_t Perturbation{0.2e-3};
 };
 
 class Swarm {
 public:
-    Swarm(CoordinateSpace ss, std::function<std::map<std::string,double_t>(Particle)> go, size_t Np, double_t otol, size_t maxit) : Space{ss}, GlobalObjective{go}, ObjectiveTolerance{otol}, MaximumIterations{maxit} {
+    Swarm(CoordinateSpace ss, std::function<ObjectiveMap(Particle)> go, size_t Np, double_t otol, size_t maxit) : Space{ss}, GlobalObjective{go}, ObjectiveTolerance{otol}, MaximumIterations{maxit} {
         std::random_device rd;
         RNG = std::mt19937{rd()};
 
@@ -139,6 +183,8 @@ public:
 
     void run();
 
+    void mpirun();
+
     std::vector<Particle> Particles;
 
     double_t Omega{1.0 / (1.0 + M_SQRT2)};
@@ -152,7 +198,7 @@ public:
 protected:
     CoordinateSpace Space;
 
-    std::function<std::map<std::string,double_t>(Particle)> GlobalObjective;
+    std::function<ObjectiveMap(Particle)> GlobalObjective;
 
     std::mt19937 RNG;
 };

test/CMakeLists.txt

@@ -46,8 +46,34 @@ set(SOURCE_FILES
 
         Model\ Templates/test_ModelTemplates.hpp
         Model\ Templates/test_DistributedWindingStator.cpp
-        Model\ Templates/test_SynchronousReluctanceRotor.cpp "Model Templates/test_ModelTemplates.hpp" "Model Templates/test_SynchronousReluctanceMachine.cpp" Optimization/test_Optimization.hpp Optimization/test_ParticleSwarmOptimization.cpp)
+        Model\ Templates/test_SynchronousReluctanceRotor.cpp
+        Model\ Templates/test_SynchronousReluctanceMachine.cpp)
+
+set(MPI_SOURCE_FILES
+        main.cpp
+
+        Optimization/test_Optimization.hpp
+        Optimization/MPI_Particle_Swarm_Optimization.cpp
+        )
+
+set(LINK_LIBRARIES
+        gtest
+        gtest_main
+        geometry
+        sketch
+        elements
+        mesh
+        materials
+        optimization
+        quadrature
+        physics
+        model_templates
+        stdc++fs
+        ${MPI_C_LIBRARIES}
+        )
 
 add_executable(run_tests ${SOURCE_FILES})
+target_link_libraries(run_tests ${LINK_LIBRARIES})
 
-target_link_libraries(run_tests gtest gtest_main geometry sketch elements mesh materials optimization quadrature physics model_templates stdc++fs)
\ No newline at end of file
+add_executable(mpi_run_tests ${MPI_SOURCE_FILES})
+target_link_libraries(mpi_run_tests ${LINK_LIBRARIES})
\ No newline at end of file

test/Optimization/MPI_Particle_Swarm_Optimization.cpp

+#include "test_Optimization.hpp"
+
+TEST(Particle_Swarm_Optimization, mpi_swarm) {
+    CoordinateSpace cs;
+
+    cs.coordinate("x") = Coordinate{-2.0,0.0,2.0};
+    cs.coordinate("y") = Coordinate{-2.0,0.0,2.0};
+
+    std::function<ObjectiveMap(Particle p)> go = [](Particle p) {
+        ObjectiveMap output;
+
+        double_t x = p.position("x");
+        double_t y = p.position("y");
+
+        output["z0"] = (x - 0.5) * (x - 0.5) + (y - 0.5) * (y - 0.5);
+
+        return output;
+    };
+
+    Swarm swarm{cs,go,21,1e-2,50};
+
+    for (size_t i = 0; i != swarm.size(); ++i) {
+        swarm.Particles[i].local_objective() = [](ObjectiveMap go) {
+            return go["z0"];
+        };
+    }
+
+    swarm.mpirun();
+
+    int rank;
+    MPI_Comm_rank(MPI_COMM_WORLD,&rank);
+    for (size_t i = 0; i != swarm.size() && rank == 0; ++i) {
+        std::cout << swarm.Particles[i] << std::endl;
+
+        //GTEST_FATAL_FAILURE_("TODO: IMPLEMENT TESTS");
+        //GTEST_FATAL_FAILURE_("TODO: CLEAN FILE STRUCTURE");
+    }
+
+    MPI_Barrier(MPI_COMM_WORLD);
+}
+
+TEST(Particle_Swarm_Optimization, mpi_multiobjective_swarm) {
+    int rank;
+    int size;
+    MPI_Comm_rank(MPI_COMM_WORLD,&rank);
+    MPI_Comm_size(MPI_COMM_WORLD,&size);
+
+    CoordinateSpace cs;
+
+    cs.coordinate("x") = Coordinate{-2.0,0.0,2.0};
+    cs.coordinate("y") = Coordinate{-2.0,0.0,2.0};
+
+    std::function<ObjectiveMap(Particle p)> go = [](Particle p) {
+        ObjectiveMap output;
+
+        double_t x = p.position("x");
+        double_t y = p.position("y");
+
+        output["z0"] = (x - 0.5) * (x - 0.5) + (y - 0.5) * (y - 0.5);
+        output["z1"] = (x + 0.5) * (x + 0.5) + (y - 0.5) * (y - 0.5);
+
+        return output;
+    };
+
+    size_t Ns{21};
+    Swarm swarm{cs,go,Ns,1e-2,50};
+
+    for (size_t i = 0; i != swarm.size(); ++i) {
+        double_t w0 = (Ns - i - 1.0) / (Ns - 1.0);
+        double_t w1 = i / (Ns - 1.0);
+
+        swarm.Particles[i].local_objective() = [w0,w1](ObjectiveMap go) {
+            return w0 * go["z0"] + w1 * go["z1"];
+        };
+    }
+
+    swarm.mpirun();
+
+    for (size_t i = 0; i != swarm.size() && rank == 0; ++i) {
+        std::cout << swarm.Particles[i] << std::endl;
+
+        //GTEST_FATAL_FAILURE_("TODO: IMPLEMENT TESTS");
+        //GTEST_FATAL_FAILURE_("TODO: CLEAN FILE STRUCTURE");
+    }
+}
\ No newline at end of file

test/Optimization/test_ParticleSwarmOptimization.cpp

@@ -48,11 +48,11 @@ TEST(Particle_Swarm_Optimization, particle_initialization) {
 
     double_t w0{1.0/3.0};
     double_t w1{2.0/3.0};
-    p.local_objective() = [w0,w1](std::map<std::string,double_t> go) {
+    p.local_objective() = [w0,w1](ObjectiveMap go) {
         return w0*go["z0"]+w1*go["z1"];
     };
 
-    std::map<std::string,double_t> go;
+    ObjectiveMap go;
     go["z0"] = M_PI;
     go["z1"] = M_E;
 
@@ -65,8 +65,8 @@ TEST(Particle_Swarm_Optimization, single_objective_swarm) {
     cs.coordinate("x") = Coordinate{-2.0,0.0,2.0};
     cs.coordinate("y") = Coordinate{-2.0,0.0,2.0};
 
-    std::function<std::map<std::string,double_t>(Particle p)> go = [](Particle p) {
-        std::map<std::string,double_t> output;
+    std::function<ObjectiveMap(Particle p)> go = [](Particle p) {
+        ObjectiveMap output;
 
         double_t x = p.position("x");
         double_t y = p.position("y");
@@ -79,7 +79,7 @@ TEST(Particle_Swarm_Optimization, single_objective_swarm) {
     Swarm swarm{cs,go,10,1e-2,20};
 
     for (size_t i = 0; i != swarm.size(); ++i) {
-        swarm.Particles[i].local_objective() = [](std::map<std::string,double_t> go) {
+        swarm.Particles[i].local_objective() = [](ObjectiveMap go) {
             return go["z0"];
         };
     }
@@ -101,8 +101,8 @@ TEST(Particle_Swarm_Optimization, multiobjective_swarm) {
     cs.coordinate("x") = Coordinate{-2.0,0.0,2.0};
     cs.coordinate("y") = Coordinate{-2.0,0.0,2.0};
 
-    std::function<std::map<std::string,double_t>(Particle p)> go = [](Particle p) {
-        std::map<std::string,double_t> output;
+    std::function<ObjectiveMap(Particle p)> go = [](Particle p) {
+        ObjectiveMap output;
 
         double_t x = p.position("x");
         double_t y = p.position("y");
@@ -119,7 +119,7 @@ TEST(Particle_Swarm_Optimization, multiobjective_swarm) {
         double_t w0 = (10 - i) / 10.0;
         double_t w1 = i / 10.0;
 
-        swarm.Particles[i].local_objective() = [w0,w1](std::map<std::string,double_t> go) {
+        swarm.Particles[i].local_objective() = [w0,w1](ObjectiveMap go) {
             return w0 * go["z0"] + w1 * go["z1"];
         };
     }