sccomp/source/solver/DifferentialEvolution.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 #ifndef INCLUDED_SCCOMP_SOURCE_DIFFERENTIALEVOLUTION_HXX
  12 #define INCLUDED_SCCOMP_SOURCE_DIFFERENTIALEVOLUTION_HXX
  13
  14 #include <vector>
  15 #include <random>
  16 #include <limits>
  17
  18 struct Individual
  19 {
  20     std::vector<double> mVariables;
  21 };
  22
  23 template <typename DataProvider> class DifferentialEvolutionAlgorithm
  24 {
  25     static constexpr double mnDifferentialWeight = 0.5; // [0, 2]
  26     static constexpr double mnCrossoverProbability = 0.9; // [0, 1]
  27
  28     static constexpr double constAcceptedPrecision = 0.000000001;
  29
  30     DataProvider& mrDataProvider;
  31
  32     size_t const mnPopulationSize;
  33     std::vector<Individual> maPopulation;
  34
  35     std::random_device maRandomDevice;
  36     std::mt19937 maGenerator;
  37     size_t const mnDimensionality;
  38
  39     std::uniform_int_distribution<> maRandomPopulation;
  40     std::uniform_int_distribution<> maRandomDimensionality;
  41     std::uniform_real_distribution<> maRandom01;
  42
  43     Individual maBestCandidate;
  44     double mfBestFitness;
  45     int mnGeneration;
  46     int mnLastChange;
  47
  48 public:
  49     DifferentialEvolutionAlgorithm(DataProvider& rDataProvider, size_t nPopulationSize)
  50         : mrDataProvider(rDataProvider)
  51         , mnPopulationSize(nPopulationSize)
  52         , maGenerator(maRandomDevice())
  53         , mnDimensionality(mrDataProvider.getDimensionality())
  54         , maRandomPopulation(0, mnPopulationSize - 1)
  55         , maRandomDimensionality(0, mnDimensionality - 1)
  56         , maRandom01(0.0, 1.0)
  57         , mfBestFitness(std::numeric_limits<double>::lowest())
  58         , mnGeneration(0)
  59         , mnLastChange(0)
  60     {
  61     }
  62
  63     std::vector<double> const& getResult() { return maBestCandidate.mVariables; }
  64
  65     int getGeneration() { return mnGeneration; }
  66
  67     int getLastChange() { return mnLastChange; }
  68
  69     void initialize()
  70     {
  71         mnGeneration = 0;
  72         mnLastChange = 0;
  73         maPopulation.clear();
  74         maBestCandidate.mVariables.clear();
  75
  76         // Initialize population with individuals that have been initialized with uniform random
  77         // noise
  78         // uniform noise means random value inside your search space
  79         maPopulation.reserve(mnPopulationSize);
  80         for (size_t i = 0; i < mnPopulationSize; ++i)
  81         {
  82             maPopulation.emplace_back();
  83             Individual& rIndividual = maPopulation.back();
  84             mrDataProvider.initializeVariables(rIndividual.mVariables, maGenerator);
  85         }
  86     }
  87
  88     // Calculate one generation
  89     bool next()
  90     {
  91         bool bBestChanged = false;
  92
  93         for (size_t agentIndex = 0; agentIndex < mnPopulationSize; ++agentIndex)
  94         {
  95             // calculate new candidate solution
  96
  97             // pick random point from population
  98             size_t x = agentIndex; // randomPopulation(generator);
  99             size_t a, b, c;
 100
 101             // create a copy of chosen random agent in population
 102             Individual& rOriginal = maPopulation[x];
 103             Individual aCandidate(rOriginal);
 104
 105             // pick three different random points from population
 106             do
 107             {
 108                 a = maRandomPopulation(maGenerator);
 109             } while (a == x);
 110
 111             do
 112             {
 113                 b = maRandomPopulation(maGenerator);
 114             } while (b == x || b == a);
 115
 116             do
 117             {
 118                 c = maRandomPopulation(maGenerator);
 119
 120             } while (c == x || c == a || c == b);
 121
 122             size_t randomIndex = maRandomDimensionality(maGenerator);
 123
 124             for (size_t index = 0; index < mnDimensionality; ++index)
 125             {
 126                 double randomCrossoverProbability = maRandom01(maGenerator);
 127                 if (index == randomIndex || randomCrossoverProbability < mnCrossoverProbability)
 128                 {
 129                     double fVarA = maPopulation[a].mVariables[index];
 130                     double fVarB = maPopulation[b].mVariables[index];
 131                     double fVarC = maPopulation[c].mVariables[index];
 132
 133                     double fNewValue = fVarA + mnDifferentialWeight * (fVarB - fVarC);
 134                     fNewValue = mrDataProvider.boundVariable(index, fNewValue);
 135                     aCandidate.mVariables[index] = fNewValue;
 136                 }
 137             }
 138
 139             double fCandidateFitness = mrDataProvider.calculateFitness(aCandidate.mVariables);
 140
 141             // see if is better than original, if so replace
 142             if (fCandidateFitness > mrDataProvider.calculateFitness(rOriginal.mVariables))
 143             {
 144                 maPopulation[x] = aCandidate;
 145
 146                 if (fCandidateFitness > mfBestFitness)
 147                 {
 148                     if (std::abs(fCandidateFitness - mfBestFitness) > constAcceptedPrecision)
 149                     {
 150                         bBestChanged = true;
 151                         mnLastChange = mnGeneration;
 152                     }
 153                     mfBestFitness = fCandidateFitness;
 154                     maBestCandidate = maPopulation[x];
 155                 }
 156             }
 157         }
 158         mnGeneration++;
 159         return bBestChanged;
 160     }
 161 };
 162
 163 #endif
 164
 165 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */