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
21#include <ql/math/optimization/differentialevolution.hpp>
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
47 EndCriteria::Type DifferentialEvolution::minimize(Problem& p, const EndCriteria& endCriteria) {
48 EndCriteria::Type ecType;
49 p.reset();
50
51 if (configuration().upperBound.empty()) {
52 upperBound_ = p.constraint().upperBound(params: p.currentValue());
53 } else {
54 QL_REQUIRE(configuration().upperBound.size() == p.currentValue().size(),
55 "wrong upper bound size in differential evolution configuration");
56 upperBound_ = configuration().upperBound;
57 }
58 if (configuration().lowerBound.empty()) {
59 lowerBound_ = p.constraint().lowerBound(params: p.currentValue());
60 } else {
61 QL_REQUIRE(configuration().lowerBound.size() == p.currentValue().size(),
62 "wrong lower bound size in differential evolution configuration");
63 lowerBound_ = configuration().lowerBound;
64 }
65 currGenSizeWeights_ =
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(new_size: 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(x: population[i].values);
78 }
79 } else {
80 population = std::vector<Candidate>(configuration().populationMembers,
81 Candidate(p.currentValue().size()));
82 fillInitialPopulation(population, p);
83 }
84
85 std::partial_sort(first: population.begin(), middle: population.begin() + 1, last: population.end(),
86 comp: 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: iteration++, ecType)) {
93 calculateNextGeneration(population, costFunction&: p);
94 std::partial_sort(first: population.begin(), middle: population.begin() + 1, last: population.end(),
95 comp: 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, statStateIterations&: stationaryPointIteration,
100 ecType))
101 break;
102 fxOld = fxNew;
103 };
104 p.setCurrentValue(bestMemberEver_.values);
105 p.setFunctionValue(bestMemberEver_.cost);
106 return ecType;
107 }
108
109 void DifferentialEvolution::calculateNextGeneration(
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(begin: population.begin(), end: population.end(), rng: rng_);
120 std::vector<Candidate> shuffledPop1 = population;
121 randomize(begin: population.begin(), end: population.end(), rng: rng_);
122 std::vector<Candidate> shuffledPop2 = population;
123 randomize(begin: population.begin(), end: population.end(), rng: rng_);
124 mirrorPopulation = shuffledPop1;
125
126 for (Size popIter = 0; popIter < population.size(); popIter++) {
127 population[popIter].values = population[popIter].values
128 + configuration().stepsizeWeight
129 * (shuffledPop1[popIter].values - shuffledPop2[popIter].values);
130 }
131 }
132 break;
133
134 case BestMemberWithJitter: {
135 randomize(begin: population.begin(), end: population.end(), rng: rng_);
136 std::vector<Candidate> shuffledPop1 = population;
137 randomize(begin: population.begin(), end: population.end(), rng: 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(),
149 bestMemberEver_);
150 }
151 break;
152
153 case CurrentToBest2Diffs: {
154 randomize(begin: population.begin(), end: population.end(), rng: rng_);
155 std::vector<Candidate> shuffledPop1 = population;
156 randomize(begin: population.begin(), end: population.end(), rng: rng_);
157
158 for (Size popIter = 0; popIter < population.size(); popIter++) {
159 population[popIter].values = oldPopulation[popIter].values
160 + configuration().stepsizeWeight
161 * (bestMemberEver_.values - oldPopulation[popIter].values)
162 + configuration().stepsizeWeight
163 * (population[popIter].values - shuffledPop1[popIter].values);
164 }
165 mirrorPopulation = shuffledPop1;
166 }
167 break;
168
169 case Rand1DiffWithPerVectorDither: {
170 randomize(begin: population.begin(), end: population.end(), rng: rng_);
171 std::vector<Candidate> shuffledPop1 = population;
172 randomize(begin: population.begin(), end: population.end(), rng: rng_);
173 std::vector<Candidate> shuffledPop2 = population;
174 randomize(begin: population.begin(), end: population.end(), rng: 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() +
179 configuration().stepsizeWeight;
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(begin: population.begin(), end: population.end(), rng: rng_);
189 std::vector<Candidate> shuffledPop1 = population;
190 randomize(begin: population.begin(), end: population.end(), rng: rng_);
191 std::vector<Candidate> shuffledPop2 = population;
192 randomize(begin: population.begin(), end: population.end(), rng: rng_);
193 mirrorPopulation = shuffledPop1;
194 Real FWeight = (1.0 - configuration().stepsizeWeight) * rng_.nextReal()
195 + configuration().stepsizeWeight;
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
203 case EitherOrWithOptimalRecombination: {
204 randomize(begin: population.begin(), end: population.end(), rng: rng_);
205 std::vector<Candidate> shuffledPop1 = population;
206 randomize(begin: population.begin(), end: population.end(), rng: rng_);
207 std::vector<Candidate> shuffledPop2 = population;
208 randomize(begin: population.begin(), end: population.end(), rng: 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
214 + configuration().stepsizeWeight
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
229 case Rand1SelfadaptiveWithRotation: {
230 randomize(begin: population.begin(), end: population.end(), rng: rng_);
231 std::vector<Candidate> shuffledPop1 = population;
232 randomize(begin: population.begin(), end: population.end(), rng: rng_);
233 std::vector<Candidate> shuffledPop2 = population;
234 randomize(begin: population.begin(), end: population.end(), rng: rng_);
235 mirrorPopulation = shuffledPop1;
236
237 adaptSizeWeights();
238
239 for (Size popIter = 0; popIter < population.size(); popIter++) {
240 if (rng_.nextReal() < 0.1){
241 population[popIter].values = rotateArray(inputArray: 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, mutantPopulation: population, mirrorPopulation, costFunction&: p);
257 }
258
259 void DifferentialEvolution::crossover(
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) {
267 adaptCrossover();
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(x: population[popIter].values);
299 } catch (Error&) {
300 population[popIter].cost = QL_MAX_REAL;
301 }
302 if (!std::isfinite(x: population[popIter].cost))
303 population[popIter].cost = QL_MAX_REAL;
304
305 }
306 }
307
308 void DifferentialEvolution::getCrossoverMask(
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
323 Array DifferentialEvolution::getMutationProbabilities(
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(x: currGenCrossover_[coIter],
338 y: (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
351 Array DifferentialEvolution::rotateArray(Array a) const {
352 randomize(begin: a.begin(), end: a.end(), rng: rng_);
353 return a;
354 }
355
356 void DifferentialEvolution::adaptSizeWeights() const {
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
370 void DifferentialEvolution::adaptCrossover() const {
371 Real crossoverChangeProb = 0.1; // [=tau2]
372 for (Real& coIter : currGenCrossover_) {
373 if (rng_.nextReal() < crossoverChangeProb)
374 coIter = rng_.nextReal();
375 }
376 }
377
378 void DifferentialEvolution::fillInitialPopulation(
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(x: 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(x: population[j].values);
392 if (!std::isfinite(x: population[j].cost))
393 population[j].cost = QL_MAX_REAL;
394 }
395 }
396
397}
398

source code of quantlib/ql/math/optimization/differentialevolution.cpp