From 0522f8b993cffe1bd56033fc1df27651805d5d80 Mon Sep 17 00:00:00 2001 From: Lars Ivar Hatledal Date: Fri, 12 Apr 2024 14:37:29 +0200 Subject: [PATCH] add PSO --- examples/libs/CMakeLists.txt | 1 + .../optimization/DifferentialEvolution.hpp | 14 ++- examples/libs/optimization/Optimizer.hpp | 5 +- .../ParticleSwarmOptimization.hpp | 114 ++++++++++++++++++ examples/libs/optimization/Problem.hpp | 43 +++++-- examples/misc/lut.cpp | 25 ++-- 6 files changed, 172 insertions(+), 30 deletions(-) create mode 100644 examples/libs/optimization/ParticleSwarmOptimization.hpp diff --git a/examples/libs/CMakeLists.txt b/examples/libs/CMakeLists.txt index c2a2a311..b827e01b 100644 --- a/examples/libs/CMakeLists.txt +++ b/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" ) diff --git a/examples/libs/optimization/DifferentialEvolution.hpp b/examples/libs/optimization/DifferentialEvolution.hpp index beb5cda4..333338bc 100644 --- a/examples/libs/optimization/DifferentialEvolution.hpp +++ b/examples/libs/optimization/DifferentialEvolution.hpp @@ -9,6 +9,14 @@ 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++) { @@ -16,10 +24,6 @@ class DifferentialEvolution: public Optimizer { } } - [[nodiscard]] const std::vector& getPopulation() const { - return population; - } - const Candidate& step(const Problem& problem) override { std::pair best{0, std::numeric_limits::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()}; diff --git a/examples/libs/optimization/Optimizer.hpp b/examples/libs/optimization/Optimizer.hpp index 6a6fec7f..09eaabc8 100644 --- a/examples/libs/optimization/Optimizer.hpp +++ b/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 diff --git a/examples/libs/optimization/ParticleSwarmOptimization.hpp b/examples/libs/optimization/ParticleSwarmOptimization.hpp new file mode 100644 index 00000000..4929b7b5 --- /dev/null +++ b/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& 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 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 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 particles; + + [[nodiscard]] const Particle& getBest() const { + std::pair best{0, std::numeric_limits::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 diff --git a/examples/libs/optimization/Problem.hpp b/examples/libs/optimization/Problem.hpp index 90e59b09..474428d8 100644 --- a/examples/libs/optimization/Problem.hpp +++ b/examples/libs/optimization/Problem.hpp @@ -8,21 +8,30 @@ #include #include -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& 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(dimensions, bounds)) {} + Problem(int dimensions, const Constraint& bounds) + : Problem(dimensions, std::vector(dimensions, bounds)) {} - Problem(int dimensions, std::vector bounds) + Problem(int dimensions, std::vector 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> solutions() const = 0; @@ -134,7 +151,7 @@ class Problem { private: int dimensions_; - std::vector bounds_; + std::vector bounds_; }; #endif//THREEPP_PROBLEM_HPP diff --git a/examples/misc/lut.cpp b/examples/misc/lut.cpp index c4b32363..0717f267 100644 --- a/examples/misc/lut.cpp +++ b/examples/misc/lut.cpp @@ -8,6 +8,7 @@ #include #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(50, 0.2f, 0.9f); + auto algorithm = std::make_unique(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)); }