add PSO

examples/libs/CMakeLists.txt

@@ -31,6 +31,7 @@ target_include_directories(pathfinding INTERFACE "${CMAKE_CURRENT_SOURCE_DIR}")
 
 add_library(optimization INTERFACE
         "${CMAKE_CURRENT_SOURCE_DIR}/optimization/DifferentialEvolution.hpp"
+        "${CMAKE_CURRENT_SOURCE_DIR}/optimization/ParticleSwarmOptimization.hpp"
         "${CMAKE_CURRENT_SOURCE_DIR}/optimization/Problem.hpp"
         "${CMAKE_CURRENT_SOURCE_DIR}/optimization/Optimizer.hpp"
 )

examples/libs/optimization/DifferentialEvolution.hpp

@@ -9,17 +9,21 @@ class DifferentialEvolution: public Optimizer {
     explicit DifferentialEvolution(size_t pop_size = 100, float mut_factor = 0.2, float cross_prob = 0.9)
         : populationSize(pop_size), differentialWeight(mut_factor), crossoverProbability(cross_prob) {}
 
+    [[nodiscard]] size_t size() const override {
+        return populationSize;
+    }
+
+    [[nodiscard]] const Candidate& getCandiateAt(size_t index) const override {
+        return population.at(index);
+    }
+
     void init(const Problem& problem) override {
         population.clear();
         for (int i = 0; i < populationSize; i++) {
             population.emplace_back(problem.generateCandidate());
         }
     }
 
-    [[nodiscard]] const std::vector<Candidate>& getPopulation() const {
-        return population;
-    }
-
     const Candidate& step(const Problem& problem) override {
 
         std::pair<int, float> best{0, std::numeric_limits<float>::max()};
@@ -28,7 +32,7 @@ class DifferentialEvolution: public Optimizer {
             const auto mutatedData = mutate(i);
             const auto mutatedCost = problem.eval(mutatedData);
             if (mutatedCost < candidate.cost()) {
-                candidate = {mutatedData, mutatedCost};
+                candidate.set(mutatedData, mutatedCost);
             }
             if (candidate.cost() < best.second) {
                 best = {i, candidate.cost()};

examples/libs/optimization/Optimizer.hpp

@@ -8,12 +8,15 @@ class Optimizer {
 
 public:
 
+    [[nodiscard]] virtual size_t size() const = 0;
+
+    [[nodiscard]] virtual const Candidate& getCandiateAt(size_t index) const = 0;
+
     virtual void init(const Problem& problem) = 0;
 
     virtual const Candidate& step(const Problem& problem) = 0;
 
     virtual ~Optimizer() = default;
-
 };
 
 #endif//THREEPP_OPTIMIMIZER_HPP

examples/libs/optimization/ParticleSwarmOptimization.hpp

@@ -0,0 +1,114 @@
+
+#ifndef THREEPP_PARTICLESWARMOPTIMIZATION_HPP
+#define THREEPP_PARTICLESWARMOPTIMIZATION_HPP
+
+#include "Optimizer.hpp"
+
+class Particle {
+public:
+    Candidate position;
+    Candidate localBestPosition;
+
+    explicit Particle(const Candidate& candidate)
+        : position(candidate), localBestPosition(candidate), velocity(candidate.size()) {
+
+        for (auto& v : velocity) {
+            v = rand() * 2 - 1;
+        }
+    }
+
+    void update(float omega, float c1, float c2, const std::vector<float>& globalBest, const Problem& problem) {
+
+        for (int i = 0; i < position.size(); i++) {
+            float r1 = rand();
+            float r2 = rand();
+
+            velocity[i] = omega * velocity[i] + c1 * r1 * (localBestPosition[i] - position[i]) + c2 * r2 * (globalBest[i] - position[i]);
+            velocity[i] = std::clamp(velocity[i], velocityBounds_.lower(), velocityBounds_.upper());
+
+            position[i] += velocity[i];
+            position[i] = std::clamp(position[i], 0.0f, 1.0f);
+        }
+
+        float cost = problem.eval(position.data());
+        position.setCost(cost);
+        if (cost < localBestPosition.cost()) {
+            localBestPosition = position;
+        }
+    }
+
+private:
+    std::vector<float> velocity;
+    Constraint velocityBounds_{-0.5, 0.5};
+
+    static float rand() {
+        static std::random_device rd;
+        static std::mt19937 gen(rd());
+        static std::uniform_real_distribution<float> randDis{0, 1};
+
+        return randDis(gen);
+    }
+};
+
+class ParticleSwarmOptimization: public Optimizer {
+
+public:
+    explicit ParticleSwarmOptimization(size_t populationSize = 30, float omega = 0.9, float c1 = 1.41, float c2 = 1.41)
+        : populationSize_(populationSize), omega_(omega), c1_(c1), c2_(c2) {}
+
+    [[nodiscard]] size_t size() const override {
+
+        return populationSize_;
+    }
+
+    [[nodiscard]] const Candidate& getCandiateAt(size_t index) const override {
+
+        return particles.at(index).localBestPosition;
+    }
+
+    void init(const Problem& problem) override {
+        particles.clear();
+
+        for (int i = 0; i < populationSize_; i++) {
+
+            particles.emplace_back(problem.generateCandidate());
+        }
+    }
+
+    const Candidate& step(const Problem& problem) override {
+
+        auto* best = &getBest();
+
+        for (int i = 0; i < populationSize_; i++) {
+            Particle& particle = particles[i];
+
+            particle.update(omega_, c1_, c2_, best->position.data(), problem);
+            if (particle.localBestPosition.cost() < best->localBestPosition.cost()) {
+                best = &particle;
+            }
+        }
+
+        return best->localBestPosition;
+    }
+
+private:
+    size_t populationSize_;
+    float omega_;
+    float c1_;
+    float c2_;
+
+    std::vector<Particle> particles;
+
+    [[nodiscard]] const Particle& getBest() const {
+        std::pair<int, float> best{0, std::numeric_limits<float>::max()};
+        for (int i = 0; i < populationSize_; i++) {
+            if (particles[i].localBestPosition.cost() < best.second) {
+                best = {i, particles[i].localBestPosition.cost()};
+            }
+        }
+
+        return particles[best.first];
+    }
+};
+
+#endif//THREEPP_PARTICLESWARMOPTIMIZATION_HPP

examples/libs/optimization/Problem.hpp

@@ -8,21 +8,30 @@
 #include <cmath>
 #include <random>
 
-class Range {
-    float lower;
-    float upper;
+class Constraint {
 
 public:
-    Range(float lower, float upper)
-        : lower(lower), upper(upper) {}
+    Constraint(float lower, float upper)
+        : lower_(lower), upper_(upper) {}
+
+    [[nodiscard]] float lower() const {
+        return lower_;
+    }
+    [[nodiscard]] float upper() const {
+        return upper_;
+    }
 
     [[nodiscard]] float normalize(float value) const {
-        return (value - lower) / (upper - lower);
+        return (value - lower_) / (upper_ - lower_);
     }
 
     [[nodiscard]] float denormalize(float value) const {
-        return (value) * (upper - lower) / 1 + lower;
+        return (value) * (upper_ - lower_) / 1 + lower_;
     }
+
+private:
+    float lower_;
+    float upper_;
 };
 
 class Candidate {
@@ -55,6 +64,15 @@ class Candidate {
         return candidate_;
     }
 
+    void set(const std::vector<float>& candidate, float cost) {
+        candidate_ = candidate;
+        cost_ = cost;
+    }
+
+    void setCost(float cost) {
+        cost_ = cost;
+    }
+
     bool operator<(const Candidate& other) const {
 
         return cost_ < other.cost_;
@@ -68,10 +86,10 @@ class Candidate {
 class Problem {
 
 public:
-    Problem(int dimensions, const Range& bounds)
-        : Problem(dimensions, std::vector<Range>(dimensions, bounds)) {}
+    Problem(int dimensions, const Constraint& bounds)
+        : Problem(dimensions, std::vector<Constraint>(dimensions, bounds)) {}
 
-    Problem(int dimensions, std::vector<Range> bounds)
+    Problem(int dimensions, std::vector<Constraint> bounds)
         : dimensions_(dimensions),
           bounds_(std::move(bounds)) {}
 
@@ -120,9 +138,8 @@ class Problem {
         for (auto& value : candidate) {
             value = dis(gen);
         }
-        auto fitness = eval(candidate);
 
-        return {candidate, fitness};
+        return {candidate,  eval(candidate)};
     }
 
     [[nodiscard]] virtual std::vector<std::vector<float>> solutions() const = 0;
@@ -134,7 +151,7 @@ class Problem {
 
 private:
     int dimensions_;
-    std::vector<Range> bounds_;
+    std::vector<Constraint> bounds_;
 };
 
 #endif//THREEPP_PROBLEM_HPP

examples/misc/lut.cpp

@@ -8,6 +8,7 @@
 #include <iostream>
 
 #include "optimization/DifferentialEvolution.hpp"
+#include "optimization/ParticleSwarmOptimization.hpp"
 
 using namespace threepp;
 
@@ -257,7 +258,8 @@ namespace {
 
 int main() {
 
-    DifferentialEvolution de(50, 0.2, 0.9);
+//    auto algorithm = std::make_unique<DifferentialEvolution>(50, 0.2f, 0.9f);
+    auto algorithm = std::make_unique<ParticleSwarmOptimization>(30, 0.9f, 1.41f, 1.41f);
 
     Lut::addColorMap("rainbow", {{0.f, 0x0000ff}, {0.001f, 0x00ffff}, {0.02f, 0xffff00}, {0.2f, 0xff0000}, {1.f, Color::darkred}});
 
@@ -317,7 +319,7 @@ int main() {
         minmax = normalizeAndApplyLut(*planeGeometry);
         normalizeAndApplyLut(*planeGeometry2);
 
-        de.init(func);
+        algorithm->init(func);
         auto solutions = func.solutions();
         solutionMesh->setCount(solutions.size());
         for (int i = 0; i < solutions.size(); i++) {
@@ -334,15 +336,15 @@ int main() {
 
     auto instancedMesh = InstancedMesh::create(
             SphereGeometry::create(0.05),
-            MeshBasicMaterial::create({{"color", Color::black}}), de.getPopulation().size());
+            MeshBasicMaterial::create({{"color", Color::black}}), algorithm->size());
     scene.add(instancedMesh);
 
     float searchSpeed = 0.7;
     ImguiFunctionalContext ui(canvas.windowPtr(), [&] {
         ImGui::SetNextWindowPos({0, 0}, 0, {0, 0});
-        ImGui::SetNextWindowSize({230, 0}, 0);
+        ImGui::SetNextWindowSize({260, 0}, 0);
 
-        ImGui::Begin("Lut");
+        ImGui::Begin("Optimization");
         ImGui::SliderFloat("Search speed", &searchSpeed, 0.1, 1);
         if (ImGui::BeginCombo("Functions", selectedFunction.c_str())) {
             for (const auto& [name, functor] : functions) {
@@ -415,12 +417,13 @@ int main() {
         ui.render();
 
         if (clock.getElapsedTime() > (1-searchSpeed)) {
-            de.step(*functions[selectedFunction]);
-            const auto& population = de.getPopulation();
-            for (auto i = 0; i < population.size(); i++) {
-                float x = math::mapLinear(population[i].data()[0], 0, 1, -gridSize / 2, gridSize / 2);
-                float z = math::mapLinear(population[i].data()[1], 0, 1, -gridSize / 2, gridSize / 2);
-                float y = math::mapLinear(functions[selectedFunction]->eval({population[i].data()[0], population[i].data()[1]}), minmax.first, minmax.second, 0, maxHeight);
+            algorithm->step(*functions[selectedFunction]);
+
+            for (auto i = 0; i < algorithm->size(); i++) {
+                const auto& candidate = algorithm->getCandiateAt(i).data();
+                float x = math::mapLinear(candidate[0], 0, 1, -gridSize / 2, gridSize / 2);
+                float z = math::mapLinear(candidate[1], 0, 1, -gridSize / 2, gridSize / 2);
+                float y = math::mapLinear(functions[selectedFunction]->eval({candidate[0], candidate[1]}), minmax.first, minmax.second, 0, maxHeight);
 
                 instancedMesh->setMatrixAt(i, m.setPosition(x, y, z));
             }