sccomp/source/solver/ParticelSwarmOptimization.hxx

   1 /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
   2 /*
   3  * This file is part of the LibreOffice project.
   4  *
   5  * This Source Code Form is subject to the terms of the Mozilla Public
   6  * License, v. 2.0. If a copy of the MPL was not distributed with this
   7  * file, You can obtain one at http://mozilla.org/MPL/2.0/.
   8  *
   9  */
  10
  11 #pragma once
  12
  13 #include <vector>
  14 #include <random>
  15 #include <limits>
  16
  17 struct Particle
  18 {
  19     Particle(size_t nDimensionality)
  20         : mVelocity(nDimensionality)
  21         , mPosition(nDimensionality)
  22         , mCurrentFitness(std::numeric_limits<float>::lowest())
  23         , mBestPosition(nDimensionality)
  24         , mBestFitness(std::numeric_limits<float>::lowest())
  25     {
  26     }
  27
  28     std::vector<double> mVelocity;
  29
  30     std::vector<double> mPosition;
  31     double mCurrentFitness;
  32
  33     std::vector<double> mBestPosition;
  34     double mBestFitness;
  35 };
  36
  37 template <typename DataProvider> class ParticleSwarmOptimizationAlgorithm
  38 {
  39 private:
  40     // inertia
  41     static constexpr double constWeight = 0.729;
  42     // cognitive coefficient
  43     static constexpr double c1 = 1.49445;
  44     // social  coefficient
  45     static constexpr double c2 = 1.49445;
  46
  47     static constexpr double constAcceptedPrecision = 0.000000001;
  48
  49     DataProvider& mrDataProvider;
  50
  51     size_t mnNumOfParticles;
  52
  53     std::vector<Particle> maSwarm;
  54
  55     std::random_device maRandomDevice;
  56     std::mt19937 maGenerator;
  57     size_t mnDimensionality;
  58
  59     std::uniform_real_distribution<> maRandom01;
  60
  61     std::vector<double> maBestPosition;
  62     double mfBestFitness;
  63     int mnGeneration;
  64     int mnLastChange;
  65
  66 public:
  67     ParticleSwarmOptimizationAlgorithm(DataProvider& rDataProvider, size_t nNumOfParticles)
  68         : mrDataProvider(rDataProvider)
  69         , mnNumOfParticles(nNumOfParticles)
  70         , maGenerator(maRandomDevice())
  71         , mnDimensionality(mrDataProvider.getDimensionality())
  72         , maRandom01(0.0, 1.0)
  73         , maBestPosition(mnDimensionality)
  74         , mfBestFitness(std::numeric_limits<float>::lowest())
  75         , mnGeneration(0)
  76         , mnLastChange(0)
  77     {
  78     }
  79
  80     std::vector<double> const& getResult() { return maBestPosition; }
  81
  82     int getGeneration() { return mnGeneration; }
  83
  84     int getLastChange() { return mnLastChange; }
  85
  86     void initialize()
  87     {
  88         mnGeneration = 0;
  89         mnLastChange = 0;
  90         maSwarm.clear();
  91
  92         mfBestFitness = std::numeric_limits<float>::lowest();
  93
  94         maSwarm.reserve(mnNumOfParticles);
  95         for (size_t i = 0; i < mnNumOfParticles; i++)
  96         {
  97             maSwarm.emplace_back(mnDimensionality);
  98             Particle& rParticle = maSwarm.back();
  99
 100             mrDataProvider.initializeVariables(rParticle.mPosition, maGenerator);
 101             mrDataProvider.initializeVariables(rParticle.mVelocity, maGenerator);
 102
 103             for (size_t k = 0; k < mnDimensionality; k++)
 104             {
 105                 rParticle.mPosition[k] = mrDataProvider.clampVariable(k, rParticle.mPosition[k]);
 106             }
 107
 108             rParticle.mCurrentFitness = mrDataProvider.calculateFitness(rParticle.mPosition);
 109
 110             for (size_t k = 0; k < mnDimensionality; k++)
 111             {
 112                 rParticle.mPosition[k] = mrDataProvider.clampVariable(k, rParticle.mPosition[k]);
 113             }
 114
 115             rParticle.mBestPosition.insert(rParticle.mBestPosition.begin(),
 116                                            rParticle.mPosition.begin(), rParticle.mPosition.end());
 117             rParticle.mBestFitness = rParticle.mCurrentFitness;
 118
 119             if (rParticle.mCurrentFitness > mfBestFitness)
 120             {
 121                 mfBestFitness = rParticle.mCurrentFitness;
 122                 maBestPosition.insert(maBestPosition.begin(), rParticle.mPosition.begin(),
 123                                       rParticle.mPosition.end());
 124             }
 125         }
 126     }
 127
 128     bool next()
 129     {
 130         bool bBestChanged = false;
 131
 132         for (Particle& rParticle : maSwarm)
 133         {
 134             double fRandom1 = maRandom01(maGenerator);
 135             double fRandom2 = maRandom01(maGenerator);
 136
 137             for (size_t k = 0; k < mnDimensionality; k++)
 138             {
 139                 rParticle.mVelocity[k]
 140                     = (constWeight * rParticle.mVelocity[k])
 141                       + (c1 * fRandom1 * (rParticle.mBestPosition[k] - rParticle.mPosition[k]))
 142                       + (c2 * fRandom2 * (maBestPosition[k] - rParticle.mPosition[k]));
 143
 144                 mrDataProvider.clampVariable(k, rParticle.mVelocity[k]);
 145
 146                 rParticle.mPosition[k] += rParticle.mVelocity[k];
 147                 rParticle.mPosition[k] = mrDataProvider.clampVariable(k, rParticle.mPosition[k]);
 148             }
 149
 150             rParticle.mCurrentFitness = mrDataProvider.calculateFitness(rParticle.mPosition);
 151
 152             if (rParticle.mCurrentFitness > rParticle.mBestFitness)
 153             {
 154                 rParticle.mBestFitness = rParticle.mCurrentFitness;
 155                 rParticle.mBestPosition.insert(rParticle.mBestPosition.begin(),
 156                                                rParticle.mPosition.begin(),
 157                                                rParticle.mPosition.end());
 158             }
 159
 160             if (rParticle.mCurrentFitness > mfBestFitness)
 161             {
 162                 if (std::abs(rParticle.mCurrentFitness - mfBestFitness) > constAcceptedPrecision)
 163                 {
 164                     bBestChanged = true;
 165                     mnLastChange = mnGeneration;
 166                 }
 167                 maBestPosition.insert(maBestPosition.begin(), rParticle.mPosition.begin(),
 168                                       rParticle.mPosition.end());
 169                 mfBestFitness = rParticle.mCurrentFitness;
 170             }
 171         }
 172         mnGeneration++;
 173         return bBestChanged;
 174     }
 175 };
 176
 177 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */