QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
differentialevolution.cpp
Go to the documentation of this file.
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2012 Ralph Schreyer
5 Copyright (C) 2012 Mateusz Kapturski
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
22#include <algorithm>
23#include <cmath>
24
25namespace QuantLib {
26
27 namespace {
28
29 struct sort_by_cost {
30 bool operator()(const DifferentialEvolution::Candidate& left,
31 const DifferentialEvolution::Candidate& right) {
32 return left.cost < right.cost;
33 }
34 };
35
36 template <class I>
37 void randomize(I begin, I end,
38 const MersenneTwisterUniformRng& rng) {
39 Size n = static_cast<Size>(end-begin);
40 for (Size i=n-1; i>0; --i) {
41 std::swap(begin[i], begin[rng.nextInt32() % (i+1)]);
42 }
43 }
44
45 }
46
48 EndCriteria::Type ecType;
49 p.reset();
50
51 if (configuration().upperBound.empty()) {
53 } else {
54 QL_REQUIRE(configuration().upperBound.size() == p.currentValue().size(),
55 "wrong upper bound size in differential evolution configuration");
57 }
58 if (configuration().lowerBound.empty()) {
60 } else {
61 QL_REQUIRE(configuration().lowerBound.size() == p.currentValue().size(),
62 "wrong lower bound size in differential evolution configuration");
64 }
66 Array(configuration().populationMembers, configuration().stepsizeWeight);
67 currGenCrossover_ = Array(configuration().populationMembers,
68 configuration().crossoverProbability);
69
70 std::vector<Candidate> population;
71 if (!configuration().initialPopulation.empty()) {
72 population.resize(configuration().initialPopulation.size());
73 for (Size i = 0; i < population.size(); ++i) {
74 population[i].values = configuration().initialPopulation[i];
75 QL_REQUIRE(population[i].values.size() == p.currentValue().size(),
76 "wrong values size in initial population");
77 population[i].cost = p.costFunction().value(population[i].values);
78 }
79 } else {
80 population = std::vector<Candidate>(configuration().populationMembers,
82 fillInitialPopulation(population, p);
83 }
84
85 std::partial_sort(population.begin(), population.begin() + 1, population.end(),
86 sort_by_cost());
87 bestMemberEver_ = population.front();
88 Real fxOld = population.front().cost;
89 Size iteration = 0, stationaryPointIteration = 0;
90
91 // main loop - calculate consecutive emerging populations
92 while (!endCriteria.checkMaxIterations(iteration++, ecType)) {
93 calculateNextGeneration(population, p);
94 std::partial_sort(population.begin(), population.begin() + 1, population.end(),
95 sort_by_cost());
96 if (population.front().cost < bestMemberEver_.cost)
97 bestMemberEver_ = population.front();
98 Real fxNew = population.front().cost;
99 if (endCriteria.checkStationaryFunctionValue(fxOld, fxNew, stationaryPointIteration,
100 ecType))
101 break;
102 fxOld = fxNew;
103 };
106 return ecType;
107 }
108
110 std::vector<Candidate>& population,
111 Problem& p) const {
112
113 std::vector<Candidate> mirrorPopulation;
114 std::vector<Candidate> oldPopulation = population;
115
116 switch (configuration().strategy) {
117
118 case Rand1Standard: {
119 randomize(population.begin(), population.end(), rng_);
120 std::vector<Candidate> shuffledPop1 = population;
121 randomize(population.begin(), population.end(), rng_);
122 std::vector<Candidate> shuffledPop2 = population;
123 randomize(population.begin(), population.end(), rng_);
124 mirrorPopulation = shuffledPop1;
125
126 for (Size popIter = 0; popIter < population.size(); popIter++) {
127 population[popIter].values = population[popIter].values
129 * (shuffledPop1[popIter].values - shuffledPop2[popIter].values);
130 }
131 }
132 break;
133
135 randomize(population.begin(), population.end(), rng_);
136 std::vector<Candidate> shuffledPop1 = population;
137 randomize(population.begin(), population.end(), rng_);
138 Array jitter(population[0].values.size(), 0.0);
139
140 for (Size popIter = 0; popIter < population.size(); popIter++) {
141 for (Real& jitterIter : jitter) {
142 jitterIter = rng_.nextReal();
143 }
144 population[popIter].values = bestMemberEver_.values
145 + (shuffledPop1[popIter].values - population[popIter].values)
146 * (0.0001 * jitter + configuration().stepsizeWeight);
147 }
148 mirrorPopulation = std::vector<Candidate>(population.size(),
150 }
151 break;
152
153 case CurrentToBest2Diffs: {
154 randomize(population.begin(), population.end(), rng_);
155 std::vector<Candidate> shuffledPop1 = population;
156 randomize(population.begin(), population.end(), rng_);
157
158 for (Size popIter = 0; popIter < population.size(); popIter++) {
159 population[popIter].values = oldPopulation[popIter].values
161 * (bestMemberEver_.values - oldPopulation[popIter].values)
163 * (population[popIter].values - shuffledPop1[popIter].values);
164 }
165 mirrorPopulation = shuffledPop1;
166 }
167 break;
168
170 randomize(population.begin(), population.end(), rng_);
171 std::vector<Candidate> shuffledPop1 = population;
172 randomize(population.begin(), population.end(), rng_);
173 std::vector<Candidate> shuffledPop2 = population;
174 randomize(population.begin(), population.end(), rng_);
175 mirrorPopulation = shuffledPop1;
176 Array FWeight = Array(population.front().values.size(), 0.0);
177 for (Real& fwIter : FWeight)
178 fwIter = (1.0 - configuration().stepsizeWeight) * rng_.nextReal() +
180 for (Size popIter = 0; popIter < population.size(); popIter++) {
181 population[popIter].values = population[popIter].values
182 + FWeight * (shuffledPop1[popIter].values - shuffledPop2[popIter].values);
183 }
184 }
185 break;
186
187 case Rand1DiffWithDither: {
188 randomize(population.begin(), population.end(), rng_);
189 std::vector<Candidate> shuffledPop1 = population;
190 randomize(population.begin(), population.end(), rng_);
191 std::vector<Candidate> shuffledPop2 = population;
192 randomize(population.begin(), population.end(), rng_);
193 mirrorPopulation = shuffledPop1;
194 Real FWeight = (1.0 - configuration().stepsizeWeight) * rng_.nextReal()
196 for (Size popIter = 0; popIter < population.size(); popIter++) {
197 population[popIter].values = population[popIter].values
198 + FWeight * (shuffledPop1[popIter].values - shuffledPop2[popIter].values);
199 }
200 }
201 break;
202
204 randomize(population.begin(), population.end(), rng_);
205 std::vector<Candidate> shuffledPop1 = population;
206 randomize(population.begin(), population.end(), rng_);
207 std::vector<Candidate> shuffledPop2 = population;
208 randomize(population.begin(), population.end(), rng_);
209 mirrorPopulation = shuffledPop1;
210 Real probFWeight = 0.5;
211 if (rng_.nextReal() < probFWeight) {
212 for (Size popIter = 0; popIter < population.size(); popIter++) {
213 population[popIter].values = oldPopulation[popIter].values
215 * (shuffledPop1[popIter].values - shuffledPop2[popIter].values);
216 }
217 } else {
218 Real K = 0.5 * (configuration().stepsizeWeight + 1); // invariant with respect to probFWeight used
219 for (Size popIter = 0; popIter < population.size(); popIter++) {
220 population[popIter].values = oldPopulation[popIter].values
221 + K
222 * (shuffledPop1[popIter].values - shuffledPop2[popIter].values
223 - 2.0 * population[popIter].values);
224 }
225 }
226 }
227 break;
228
230 randomize(population.begin(), population.end(), rng_);
231 std::vector<Candidate> shuffledPop1 = population;
232 randomize(population.begin(), population.end(), rng_);
233 std::vector<Candidate> shuffledPop2 = population;
234 randomize(population.begin(), population.end(), rng_);
235 mirrorPopulation = shuffledPop1;
236
238
239 for (Size popIter = 0; popIter < population.size(); popIter++) {
240 if (rng_.nextReal() < 0.1){
241 population[popIter].values = rotateArray(bestMemberEver_.values);
242 }else {
243 population[popIter].values = bestMemberEver_.values
244 + currGenSizeWeights_[popIter]
245 * (shuffledPop1[popIter].values - shuffledPop2[popIter].values);
246 }
247 }
248 }
249 break;
250
251 default:
252 QL_FAIL("Unknown strategy ("
253 << Integer(configuration().strategy) << ")");
254 }
255 // in order to avoid unnecessary copying we use the same population object for mutants
256 crossover(oldPopulation, population, population, mirrorPopulation, p);
257 }
258
260 const std::vector<Candidate>& oldPopulation,
261 std::vector<Candidate>& population,
262 const std::vector<Candidate>& mutantPopulation,
263 const std::vector<Candidate>& mirrorPopulation,
264 Problem& p) const {
265
266 if (configuration().crossoverIsAdaptive) {
268 }
269
270 Array mutationProbabilities = getMutationProbabilities(population);
271
272 std::vector<Array> crossoverMask(population.size(),
273 Array(population.front().values.size(), 1.0));
274 std::vector<Array> invCrossoverMask = crossoverMask;
275 getCrossoverMask(crossoverMask, invCrossoverMask, mutationProbabilities);
276
277 // crossover of the old and mutant population
278 for (Size popIter = 0; popIter < population.size(); popIter++) {
279 population[popIter].values = oldPopulation[popIter].values * invCrossoverMask[popIter]
280 + mutantPopulation[popIter].values * crossoverMask[popIter];
281 // immediately apply bounds if specified
282 if (configuration().applyBounds) {
283 for (Size memIter = 0; memIter < population[popIter].values.size(); memIter++) {
284 if (population[popIter].values[memIter] > upperBound_[memIter])
285 population[popIter].values[memIter] = upperBound_[memIter]
286 + rng_.nextReal()
287 * (mirrorPopulation[popIter].values[memIter]
288 - upperBound_[memIter]);
289 if (population[popIter].values[memIter] < lowerBound_[memIter])
290 population[popIter].values[memIter] = lowerBound_[memIter]
291 + rng_.nextReal()
292 * (mirrorPopulation[popIter].values[memIter]
293 - lowerBound_[memIter]);
294 }
295 }
296 // evaluate objective function as soon as possible to avoid unnecessary loops
297 try {
298 population[popIter].cost = p.value(population[popIter].values);
299 } catch (Error&) {
300 population[popIter].cost = QL_MAX_REAL;
301 }
302 if (!std::isfinite(population[popIter].cost))
303 population[popIter].cost = QL_MAX_REAL;
304
305 }
306 }
307
309 std::vector<Array> & crossoverMask,
310 std::vector<Array> & invCrossoverMask,
311 const Array & mutationProbabilities) const {
312 for (Size cmIter = 0; cmIter < crossoverMask.size(); cmIter++) {
313 for (Size memIter = 0; memIter < crossoverMask[cmIter].size(); memIter++) {
314 if (rng_.nextReal() < mutationProbabilities[cmIter]) {
315 invCrossoverMask[cmIter][memIter] = 0.0;
316 } else {
317 crossoverMask[cmIter][memIter] = 0.0;
318 }
319 }
320 }
321 }
322
324 const std::vector<Candidate> & population) const {
325 Array mutationProbabilities = currGenCrossover_;
326 switch (configuration().crossoverType) {
327 case Normal:
328 break;
329 case Binomial:
330 mutationProbabilities = currGenCrossover_
331 * (1.0 - 1.0 / population.front().values.size())
332 + 1.0 / population.front().values.size();
333 break;
334 case Exponential:
335 for (Size coIter = 0;coIter< currGenCrossover_.size(); coIter++){
336 mutationProbabilities[coIter] =
337 (1.0 - std::pow(currGenCrossover_[coIter],
338 (int) population.front().values.size()))
339 / (population.front().values.size()
340 * (1.0 - currGenCrossover_[coIter]));
341 }
342 break;
343 default:
344 QL_FAIL("Unknown crossover type ("
345 << Integer(configuration().crossoverType) << ")");
346 break;
347 }
348 return mutationProbabilities;
349 }
350
352 randomize(a.begin(), a.end(), rng_);
353 return a;
354 }
355
357 // [=Fl & =Fu] respectively see Brest, J. et al., 2006,
358 // "Self-Adapting Control Parameters in Differential
359 // Evolution"
360 Real sizeWeightLowerBound = 0.1, sizeWeightUpperBound = 0.9;
361 // [=tau1] A Comparative Study on Numerical Benchmark
362 // Problems." page 649 for reference
363 Real sizeWeightChangeProb = 0.1;
364 for (Real& currGenSizeWeight : currGenSizeWeights_) {
365 if (rng_.nextReal() < sizeWeightChangeProb)
366 currGenSizeWeight = sizeWeightLowerBound + rng_.nextReal() * sizeWeightUpperBound;
367 }
368 }
369
371 Real crossoverChangeProb = 0.1; // [=tau2]
372 for (Real& coIter : currGenCrossover_) {
373 if (rng_.nextReal() < crossoverChangeProb)
374 coIter = rng_.nextReal();
375 }
376 }
377
379 std::vector<Candidate> & population,
380 const Problem& p) const {
381
382 // use initial values provided by the user
383 population.front().values = p.currentValue();
384 population.front().cost = p.costFunction().value(population.front().values);
385 // rest of the initial population is random
386 for (Size j = 1; j < population.size(); ++j) {
387 for (Size i = 0; i < p.currentValue().size(); ++i) {
388 Real l = lowerBound_[i], u = upperBound_[i];
389 population[j].values[i] = l + (u-l)*rng_.nextReal();
390 }
391 population[j].cost = p.costFunction().value(population[j].values);
392 if (!std::isfinite(population[j].cost))
393 population[j].cost = QL_MAX_REAL;
394 }
395 }
396
397}
1-D array used in linear algebra.
Definition: array.hpp:52
const_iterator end() const
Definition: array.hpp:511
Size size() const
dimension of the array
Definition: array.hpp:495
Real front() const
Definition: array.hpp:451
const_iterator begin() const
Definition: array.hpp:503
Array lowerBound(const Array &params) const
Definition: constraint.hpp:66
Array upperBound(const Array &params) const
Definition: constraint.hpp:58
virtual Real value(const Array &x) const
method to overload to compute the cost function value in x
void fillInitialPopulation(std::vector< Candidate > &population, const Problem &p) const
const Configuration & configuration() const
void getCrossoverMask(std::vector< Array > &crossoverMask, std::vector< Array > &invCrossoverMask, const Array &mutationProbabilities) const
void calculateNextGeneration(std::vector< Candidate > &population, Problem &costFunction) const
void crossover(const std::vector< Candidate > &oldPopulation, std::vector< Candidate > &population, const std::vector< Candidate > &mutantPopulation, const std::vector< Candidate > &mirrorPopulation, Problem &costFunction) const
Array getMutationProbabilities(const std::vector< Candidate > &population) const
EndCriteria::Type minimize(Problem &p, const EndCriteria &endCriteria) override
minimize the optimization problem P
Array rotateArray(Array inputArray) const
Criteria to end optimization process:
Definition: endcriteria.hpp:40
bool checkStationaryFunctionValue(Real fxOld, Real fxNew, Size &statStateIterations, EndCriteria::Type &ecType) const
Definition: endcriteria.cpp:79
bool checkMaxIterations(Size iteration, EndCriteria::Type &ecType) const
Definition: endcriteria.cpp:56
Base error class.
Definition: errors.hpp:39
Real nextReal() const
return a random number in the (0.0, 1.0)-interval
Constrained optimization problem.
Definition: problem.hpp:42
const Array & currentValue() const
current value of the local minimum
Definition: problem.hpp:81
Constraint & constraint() const
Constraint.
Definition: problem.hpp:71
Real value(const Array &x)
call cost function computation and increment evaluation counter
Definition: problem.hpp:116
void setFunctionValue(Real functionValue)
Definition: problem.hpp:83
CostFunction & costFunction() const
Cost function.
Definition: problem.hpp:74
void setCurrentValue(const Array &currentValue)
Definition: problem.hpp:76
Differential Evolution optimization method.
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
#define QL_FAIL(message)
throw an error (possibly with file and line information)
Definition: errors.hpp:92
#define QL_MAX_REAL
Definition: qldefines.hpp:176
QL_REAL Real
real number
Definition: types.hpp:50
QL_INTEGER Integer
integer number
Definition: types.hpp:35
std::size_t Size
size of a container
Definition: types.hpp:58
Definition: any.hpp:35