/src/quantlib/ql/math/optimization/differentialevolution.cpp

Source (jump to first uncovered line)
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2012 Ralph Schreyer
 Copyright (C) 2012 Mateusz Kapturski

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <http://quantlib.org/license.shtml>.

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/

#include <ql/math/optimization/differentialevolution.hpp>
#include <algorithm>
#include <cmath>

namespace QuantLib {

    namespace {

        struct sort_by_cost {
            bool operator()(const DifferentialEvolution::Candidate& left,
                            const DifferentialEvolution::Candidate& right) {
                return left.cost < right.cost;
            }
        };

        template <class I>
        void randomize(I begin, I end,
                       const MersenneTwisterUniformRng& rng) {
            Size n = static_cast<Size>(end-begin);
            for (Size i=n-1; i>0; --i) {
                std::swap(begin[i], begin[rng.nextInt32() % (i+1)]);
            }
        }

    }

    EndCriteria::Type DifferentialEvolution::minimize(Problem& p, const EndCriteria& endCriteria) {
        EndCriteria::Type ecType;
        p.reset();

        if (configuration().upperBound.empty()) {
            upperBound_ = p.constraint().upperBound(p.currentValue());
        } else {
            QL_REQUIRE(configuration().upperBound.size() == p.currentValue().size(),
                       "wrong upper bound size in differential evolution configuration");
            upperBound_ = configuration().upperBound;
        }
        if (configuration().lowerBound.empty()) {
            lowerBound_ = p.constraint().lowerBound(p.currentValue());
        } else {
            QL_REQUIRE(configuration().lowerBound.size() == p.currentValue().size(),
                       "wrong lower bound size in differential evolution configuration");
            lowerBound_ = configuration().lowerBound;
        }
        currGenSizeWeights_ =
            Array(configuration().populationMembers, configuration().stepsizeWeight);
        currGenCrossover_ = Array(configuration().populationMembers,
                                  configuration().crossoverProbability);

        std::vector<Candidate> population;
        if (!configuration().initialPopulation.empty()) {
            population.resize(configuration().initialPopulation.size());
            for (Size i = 0; i < population.size(); ++i) {
                population[i].values = configuration().initialPopulation[i];
                QL_REQUIRE(population[i].values.size() == p.currentValue().size(),
                           "wrong values size in initial population");
                population[i].cost = p.costFunction().value(population[i].values);
            }
        } else {
            population = std::vector<Candidate>(configuration().populationMembers,
                                                Candidate(p.currentValue().size()));
            fillInitialPopulation(population, p);
        }

        std::partial_sort(population.begin(), population.begin() + 1, population.end(),
                          sort_by_cost());
        bestMemberEver_ = population.front();
        Real fxOld = population.front().cost;
        Size iteration = 0, stationaryPointIteration = 0;

        // main loop - calculate consecutive emerging populations
        while (!endCriteria.checkMaxIterations(iteration++, ecType)) {
            calculateNextGeneration(population, p);
            std::partial_sort(population.begin(), population.begin() + 1, population.end(),
                              sort_by_cost());
            if (population.front().cost < bestMemberEver_.cost)
                bestMemberEver_ = population.front();
            Real fxNew = population.front().cost;
            if (endCriteria.checkStationaryFunctionValue(fxOld, fxNew, stationaryPointIteration,
                                                         ecType))
                break;
            fxOld = fxNew;
        };
        p.setCurrentValue(bestMemberEver_.values);
        p.setFunctionValue(bestMemberEver_.cost);
        return ecType;
    }

    void DifferentialEvolution::calculateNextGeneration(
                                     std::vector<Candidate>& population,
                                     Problem& p) const {

        std::vector<Candidate> mirrorPopulation;
        std::vector<Candidate> oldPopulation = population;

        switch (configuration().strategy) {

          case Rand1Standard: {
              randomize(population.begin(), population.end(), rng_);
              std::vector<Candidate> shuffledPop1 = population;
              randomize(population.begin(), population.end(), rng_);
              std::vector<Candidate> shuffledPop2 = population;
              randomize(population.begin(), population.end(), rng_);
              mirrorPopulation = shuffledPop1;

              for (Size popIter = 0; popIter < population.size(); popIter++) {
                  population[popIter].values = population[popIter].values
                      + configuration().stepsizeWeight
                      * (shuffledPop1[popIter].values - shuffledPop2[popIter].values);
              }
          }
            break;

          case BestMemberWithJitter: {
              randomize(population.begin(), population.end(), rng_);
              std::vector<Candidate> shuffledPop1 = population;
              randomize(population.begin(), population.end(), rng_);
              Array jitter(population[0].values.size(), 0.0);

              for (Size popIter = 0; popIter < population.size(); popIter++) {
                  for (Real& jitterIter : jitter) {
                      jitterIter = rng_.nextReal();
                  }
                  population[popIter].values = bestMemberEver_.values
                      + (shuffledPop1[popIter].values - population[popIter].values)
                      * (0.0001 * jitter + configuration().stepsizeWeight);
              }
              mirrorPopulation = std::vector<Candidate>(population.size(),
                                                        bestMemberEver_);
          }
            break;

          case CurrentToBest2Diffs: {
              randomize(population.begin(), population.end(), rng_);
              std::vector<Candidate> shuffledPop1 = population;
              randomize(population.begin(), population.end(), rng_);

              for (Size popIter = 0; popIter < population.size(); popIter++) {
                  population[popIter].values = oldPopulation[popIter].values
                      + configuration().stepsizeWeight
                      * (bestMemberEver_.values - oldPopulation[popIter].values)
                      + configuration().stepsizeWeight
                      * (population[popIter].values - shuffledPop1[popIter].values);
              }
              mirrorPopulation = shuffledPop1;
          }
            break;

          case Rand1DiffWithPerVectorDither: {
              randomize(population.begin(), population.end(), rng_);
              std::vector<Candidate> shuffledPop1 = population;
              randomize(population.begin(), population.end(), rng_);
              std::vector<Candidate> shuffledPop2 = population;
              randomize(population.begin(), population.end(), rng_);
              mirrorPopulation = shuffledPop1;
              Array FWeight = Array(population.front().values.size(), 0.0);
              for (Real& fwIter : FWeight)
                  fwIter = (1.0 - configuration().stepsizeWeight) * rng_.nextReal() +
                           configuration().stepsizeWeight;
              for (Size popIter = 0; popIter < population.size(); popIter++) {
                  population[popIter].values = population[popIter].values
                      + FWeight * (shuffledPop1[popIter].values - shuffledPop2[popIter].values);
              }
          }
            break;

          case Rand1DiffWithDither: {
              randomize(population.begin(), population.end(), rng_);
              std::vector<Candidate> shuffledPop1 = population;
              randomize(population.begin(), population.end(), rng_);
              std::vector<Candidate> shuffledPop2 = population;
              randomize(population.begin(), population.end(), rng_);
              mirrorPopulation = shuffledPop1;
              Real FWeight = (1.0 - configuration().stepsizeWeight) * rng_.nextReal()
                  + configuration().stepsizeWeight;
              for (Size popIter = 0; popIter < population.size(); popIter++) {
                  population[popIter].values = population[popIter].values
                      + FWeight * (shuffledPop1[popIter].values - shuffledPop2[popIter].values);
              }
          }
            break;

          case EitherOrWithOptimalRecombination: {
              randomize(population.begin(), population.end(), rng_);
              std::vector<Candidate> shuffledPop1 = population;
              randomize(population.begin(), population.end(), rng_);
              std::vector<Candidate> shuffledPop2 = population;
              randomize(population.begin(), population.end(), rng_);
              mirrorPopulation = shuffledPop1;
              Real probFWeight = 0.5;
              if (rng_.nextReal() < probFWeight) {
                  for (Size popIter = 0; popIter < population.size(); popIter++) {
                      population[popIter].values = oldPopulation[popIter].values
                          + configuration().stepsizeWeight
                          * (shuffledPop1[popIter].values - shuffledPop2[popIter].values);
                  }
              } else {
                  Real K = 0.5 * (configuration().stepsizeWeight + 1); // invariant with respect to probFWeight used
                  for (Size popIter = 0; popIter < population.size(); popIter++) {
                      population[popIter].values = oldPopulation[popIter].values
                          + K
                          * (shuffledPop1[popIter].values - shuffledPop2[popIter].values
                             - 2.0 * population[popIter].values);
                  }
              }
          }
            break;

          case Rand1SelfadaptiveWithRotation: {
              randomize(population.begin(), population.end(), rng_);
              std::vector<Candidate> shuffledPop1 = population;
              randomize(population.begin(), population.end(), rng_);
              std::vector<Candidate> shuffledPop2 = population;
              randomize(population.begin(), population.end(), rng_);
              mirrorPopulation = shuffledPop1;

              adaptSizeWeights();

              for (Size popIter = 0; popIter < population.size(); popIter++) {
                  if (rng_.nextReal() < 0.1){
                      population[popIter].values = rotateArray(bestMemberEver_.values);
                  }else {
                      population[popIter].values = bestMemberEver_.values
                          + currGenSizeWeights_[popIter]
                          * (shuffledPop1[popIter].values - shuffledPop2[popIter].values);
                  }
              }
          }
            break;

          default:
            QL_FAIL("Unknown strategy ("
                    << Integer(configuration().strategy) << ")");
        }
        // in order to avoid unnecessary copying we use the same population object for mutants
        crossover(oldPopulation, population, population, mirrorPopulation, p);
    }

    void DifferentialEvolution::crossover(
                               const std::vector<Candidate>& oldPopulation,
                               std::vector<Candidate>& population,
                               const std::vector<Candidate>& mutantPopulation,
                               const std::vector<Candidate>& mirrorPopulation,
                               Problem& p) const {

        if (configuration().crossoverIsAdaptive) {
            adaptCrossover();
        }

        Array mutationProbabilities = getMutationProbabilities(population);

        std::vector<Array> crossoverMask(population.size(),
                                         Array(population.front().values.size(), 1.0));
        std::vector<Array> invCrossoverMask = crossoverMask;
        getCrossoverMask(crossoverMask, invCrossoverMask, mutationProbabilities);

        // crossover of the old and mutant population
        for (Size popIter = 0; popIter < population.size(); popIter++) {
            population[popIter].values = oldPopulation[popIter].values * invCrossoverMask[popIter]
                + mutantPopulation[popIter].values * crossoverMask[popIter];
            // immediately apply bounds if specified
            if (configuration().applyBounds) {
                for (Size memIter = 0; memIter < population[popIter].values.size(); memIter++) {
                    if (population[popIter].values[memIter] > upperBound_[memIter])
                        population[popIter].values[memIter] = upperBound_[memIter]
                            + rng_.nextReal()
                            * (mirrorPopulation[popIter].values[memIter]
                               - upperBound_[memIter]);
                    if (population[popIter].values[memIter] < lowerBound_[memIter])
                        population[popIter].values[memIter] = lowerBound_[memIter]
                            + rng_.nextReal()
                            * (mirrorPopulation[popIter].values[memIter]
                               - lowerBound_[memIter]);
                }
            }
            // evaluate objective function as soon as possible to avoid unnecessary loops
            try {
                population[popIter].cost = p.value(population[popIter].values);
            } catch (Error&) {
                population[popIter].cost = QL_MAX_REAL;
            }
            if (!std::isfinite(population[popIter].cost))
                population[popIter].cost = QL_MAX_REAL;

        }
    }

    void DifferentialEvolution::getCrossoverMask(
                                  std::vector<Array> & crossoverMask,
                                  std::vector<Array> & invCrossoverMask,
                                  const Array & mutationProbabilities) const {
        for (Size cmIter = 0; cmIter < crossoverMask.size(); cmIter++) {
            for (Size memIter = 0; memIter < crossoverMask[cmIter].size(); memIter++) {
                if (rng_.nextReal() < mutationProbabilities[cmIter]) {
                    invCrossoverMask[cmIter][memIter] = 0.0;
                } else {
                    crossoverMask[cmIter][memIter] = 0.0;
                }
            }
        }
    }

    Array DifferentialEvolution::getMutationProbabilities(
                            const std::vector<Candidate> & population) const {
        Array mutationProbabilities = currGenCrossover_;
        switch (configuration().crossoverType) {
          case Normal:
            break;
          case Binomial:
            mutationProbabilities = currGenCrossover_
                * (1.0 - 1.0 / population.front().values.size())
                + 1.0 / population.front().values.size();
            break;
          case Exponential:
            for (Size coIter = 0;coIter< currGenCrossover_.size(); coIter++){
                mutationProbabilities[coIter] =
                    (1.0 - std::pow(currGenCrossover_[coIter],
                                    (int) population.front().values.size()))
                    / (population.front().values.size()
                       * (1.0 - currGenCrossover_[coIter]));
            }
            break;
          default:
            QL_FAIL("Unknown crossover type ("
                    << Integer(configuration().crossoverType) << ")");
            break;
        }
        return mutationProbabilities;
    }

    Array DifferentialEvolution::rotateArray(Array a) const {
        randomize(a.begin(), a.end(), rng_);
        return a;
    }

    void DifferentialEvolution::adaptSizeWeights() const {
        // [=Fl & =Fu] respectively see Brest, J. et al., 2006,
        // "Self-Adapting Control Parameters in Differential
        // Evolution"
        Real sizeWeightLowerBound = 0.1, sizeWeightUpperBound = 0.9;
         // [=tau1] A Comparative Study on Numerical Benchmark
         // Problems." page 649 for reference
        Real sizeWeightChangeProb = 0.1;
        for (Real& currGenSizeWeight : currGenSizeWeights_) {
            if (rng_.nextReal() < sizeWeightChangeProb)
                currGenSizeWeight = sizeWeightLowerBound + rng_.nextReal() * sizeWeightUpperBound;
        }
    }

    void DifferentialEvolution::adaptCrossover() const {
        Real crossoverChangeProb = 0.1; // [=tau2]
        for (Real& coIter : currGenCrossover_) {
            if (rng_.nextReal() < crossoverChangeProb)
                coIter = rng_.nextReal();
        }
    }

    void DifferentialEvolution::fillInitialPopulation(
                                          std::vector<Candidate> & population,
                                          const Problem& p) const {

        // use initial values provided by the user
        population.front().values = p.currentValue();
        population.front().cost = p.costFunction().value(population.front().values);
        // rest of the initial population is random
        for (Size j = 1; j < population.size(); ++j) {
            for (Size i = 0; i < p.currentValue().size(); ++i) {
                Real l = lowerBound_[i], u = upperBound_[i];
                population[j].values[i] = l + (u-l)*rng_.nextReal();
            }
            population[j].cost = p.costFunction().value(population[j].values);
            if (!std::isfinite(population[j].cost))
                population[j].cost = QL_MAX_REAL;
        }
    }

}

Coverage Report

Created: 2025-08-11 06:28

Line	Count	Source (jump to first uncovered line)
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
25		namespace QuantLib {
26
27		namespace {
28
29		struct sort_by_cost {
30		bool operator()(const DifferentialEvolution::Candidate& left,
31	0	const DifferentialEvolution::Candidate& right) {
32	0	return left.cost < right.cost;
33	0	}
34		};
35
36		template <class I>
37		void randomize(I begin, I end,
38	0	const MersenneTwisterUniformRng& rng) {
39	0	Size n = static_cast<Size>(end-begin);
40	0	for (Size i=n-1; i>0; --i) {
41	0	std::swap(begin[i], begin[rng.nextInt32() % (i+1)]);
42	0	}
43	0	} Unexecuted instantiation: differentialevolution.cpp:void QuantLib::(anonymous namespace)::randomize<std::__1::__wrap_iter<QuantLib::DifferentialEvolution::Candidate> >(std::__1::__wrap_iter<QuantLib::DifferentialEvolution::Candidate>, std::__1::__wrap_iter<QuantLib::DifferentialEvolution::Candidate>, QuantLib::MersenneTwisterUniformRng const&) Unexecuted instantiation: differentialevolution.cpp:void QuantLib::(anonymous namespace)::randomize<double>(double, double, QuantLib::MersenneTwisterUniformRng const&)
44
45		}
46
47	0	EndCriteria::Type DifferentialEvolution::minimize(Problem& p, const EndCriteria& endCriteria) {
48	0	EndCriteria::Type ecType;
49	0	p.reset();
50
51	0	if (configuration().upperBound.empty()) {
52	0	upperBound_ = p.constraint().upperBound(p.currentValue());
53	0	} else {
54	0	QL_REQUIRE(configuration().upperBound.size() == p.currentValue().size(),
55	0	"wrong upper bound size in differential evolution configuration");
56	0	upperBound_ = configuration().upperBound;
57	0	}
58	0	if (configuration().lowerBound.empty()) {
59	0	lowerBound_ = p.constraint().lowerBound(p.currentValue());
60	0	} else {
61	0	QL_REQUIRE(configuration().lowerBound.size() == p.currentValue().size(),
62	0	"wrong lower bound size in differential evolution configuration");
63	0	lowerBound_ = configuration().lowerBound;
64	0	}
65	0	currGenSizeWeights_ =
66	0	Array(configuration().populationMembers, configuration().stepsizeWeight);
67	0	currGenCrossover_ = Array(configuration().populationMembers,
68	0	configuration().crossoverProbability);
69
70	0	std::vector<Candidate> population;
71	0	if (!configuration().initialPopulation.empty()) {
72	0	population.resize(configuration().initialPopulation.size());
73	0	for (Size i = 0; i < population.size(); ++i) {
74	0	population[i].values = configuration().initialPopulation[i];
75	0	QL_REQUIRE(population[i].values.size() == p.currentValue().size(),
76	0	"wrong values size in initial population");
77	0	population[i].cost = p.costFunction().value(population[i].values);
78	0	}
79	0	} else {
80	0	population = std::vector<Candidate>(configuration().populationMembers,
81	0	Candidate(p.currentValue().size()));
82	0	fillInitialPopulation(population, p);
83	0	}
84
85	0	std::partial_sort(population.begin(), population.begin() + 1, population.end(),
86	0	sort_by_cost());
87	0	bestMemberEver_ = population.front();
88	0	Real fxOld = population.front().cost;
89	0	Size iteration = 0, stationaryPointIteration = 0;
90
91		// main loop - calculate consecutive emerging populations
92	0	while (!endCriteria.checkMaxIterations(iteration++, ecType)) {
93	0	calculateNextGeneration(population, p);
94	0	std::partial_sort(population.begin(), population.begin() + 1, population.end(),
95	0	sort_by_cost());
96	0	if (population.front().cost < bestMemberEver_.cost)
97	0	bestMemberEver_ = population.front();
98	0	Real fxNew = population.front().cost;
99	0	if (endCriteria.checkStationaryFunctionValue(fxOld, fxNew, stationaryPointIteration,
100	0	ecType))
101	0	break;
102	0	fxOld = fxNew;
103	0	};
104	0	p.setCurrentValue(bestMemberEver_.values);
105	0	p.setFunctionValue(bestMemberEver_.cost);
106	0	return ecType;
107	0	}
108
109		void DifferentialEvolution::calculateNextGeneration(
110		std::vector<Candidate>& population,
111	0	Problem& p) const {
112
113	0	std::vector<Candidate> mirrorPopulation;
114	0	std::vector<Candidate> oldPopulation = population;
115
116	0	switch (configuration().strategy) {
117
118	0	case Rand1Standard: {
119	0	randomize(population.begin(), population.end(), rng_);
120	0	std::vector<Candidate> shuffledPop1 = population;
121	0	randomize(population.begin(), population.end(), rng_);
122	0	std::vector<Candidate> shuffledPop2 = population;
123	0	randomize(population.begin(), population.end(), rng_);
124	0	mirrorPopulation = shuffledPop1;
125
126	0	for (Size popIter = 0; popIter < population.size(); popIter++) {
127	0	population[popIter].values = population[popIter].values
128	0	+ configuration().stepsizeWeight
129	0	* (shuffledPop1[popIter].values - shuffledPop2[popIter].values);
130	0	}
131	0	}
132	0	break;
133
134	0	case BestMemberWithJitter: {
135	0	randomize(population.begin(), population.end(), rng_);
136	0	std::vector<Candidate> shuffledPop1 = population;
137	0	randomize(population.begin(), population.end(), rng_);
138	0	Array jitter(population[0].values.size(), 0.0);
139
140	0	for (Size popIter = 0; popIter < population.size(); popIter++) {
141	0	for (Real& jitterIter : jitter) {
142	0	jitterIter = rng_.nextReal();
143	0	}
144	0	population[popIter].values = bestMemberEver_.values
145	0	+ (shuffledPop1[popIter].values - population[popIter].values)
146	0	* (0.0001 * jitter + configuration().stepsizeWeight);
147	0	}
148	0	mirrorPopulation = std::vector<Candidate>(population.size(),
149	0	bestMemberEver_);
150	0	}
151	0	break;
152
153	0	case CurrentToBest2Diffs: {
154	0	randomize(population.begin(), population.end(), rng_);
155	0	std::vector<Candidate> shuffledPop1 = population;
156	0	randomize(population.begin(), population.end(), rng_);
157
158	0	for (Size popIter = 0; popIter < population.size(); popIter++) {
159	0	population[popIter].values = oldPopulation[popIter].values
160	0	+ configuration().stepsizeWeight
161	0	* (bestMemberEver_.values - oldPopulation[popIter].values)
162	0	+ configuration().stepsizeWeight
163	0	* (population[popIter].values - shuffledPop1[popIter].values);
164	0	}
165	0	mirrorPopulation = shuffledPop1;
166	0	}
167	0	break;
168
169	0	case Rand1DiffWithPerVectorDither: {
170	0	randomize(population.begin(), population.end(), rng_);
171	0	std::vector<Candidate> shuffledPop1 = population;
172	0	randomize(population.begin(), population.end(), rng_);
173	0	std::vector<Candidate> shuffledPop2 = population;
174	0	randomize(population.begin(), population.end(), rng_);
175	0	mirrorPopulation = shuffledPop1;
176	0	Array FWeight = Array(population.front().values.size(), 0.0);
177	0	for (Real& fwIter : FWeight)
178	0	fwIter = (1.0 - configuration().stepsizeWeight) * rng_.nextReal() +
179	0	configuration().stepsizeWeight;
180	0	for (Size popIter = 0; popIter < population.size(); popIter++) {
181	0	population[popIter].values = population[popIter].values
182	0	+ FWeight * (shuffledPop1[popIter].values - shuffledPop2[popIter].values);
183	0	}
184	0	}
185	0	break;
186
187	0	case Rand1DiffWithDither: {
188	0	randomize(population.begin(), population.end(), rng_);
189	0	std::vector<Candidate> shuffledPop1 = population;
190	0	randomize(population.begin(), population.end(), rng_);
191	0	std::vector<Candidate> shuffledPop2 = population;
192	0	randomize(population.begin(), population.end(), rng_);
193	0	mirrorPopulation = shuffledPop1;
194	0	Real FWeight = (1.0 - configuration().stepsizeWeight) * rng_.nextReal()
195	0	+ configuration().stepsizeWeight;
196	0	for (Size popIter = 0; popIter < population.size(); popIter++) {
197	0	population[popIter].values = population[popIter].values
198	0	+ FWeight * (shuffledPop1[popIter].values - shuffledPop2[popIter].values);
199	0	}
200	0	}
201	0	break;
202
203	0	case EitherOrWithOptimalRecombination: {
204	0	randomize(population.begin(), population.end(), rng_);
205	0	std::vector<Candidate> shuffledPop1 = population;
206	0	randomize(population.begin(), population.end(), rng_);
207	0	std::vector<Candidate> shuffledPop2 = population;
208	0	randomize(population.begin(), population.end(), rng_);
209	0	mirrorPopulation = shuffledPop1;
210	0	Real probFWeight = 0.5;
211	0	if (rng_.nextReal() < probFWeight) {
212	0	for (Size popIter = 0; popIter < population.size(); popIter++) {
213	0	population[popIter].values = oldPopulation[popIter].values
214	0	+ configuration().stepsizeWeight
215	0	* (shuffledPop1[popIter].values - shuffledPop2[popIter].values);
216	0	}
217	0	} else {
218	0	Real K = 0.5 * (configuration().stepsizeWeight + 1); // invariant with respect to probFWeight used
219	0	for (Size popIter = 0; popIter < population.size(); popIter++) {
220	0	population[popIter].values = oldPopulation[popIter].values
221	0	+ K
222	0	* (shuffledPop1[popIter].values - shuffledPop2[popIter].values
223	0	- 2.0 * population[popIter].values);
224	0	}
225	0	}
226	0	}
227	0	break;
228
229	0	case Rand1SelfadaptiveWithRotation: {
230	0	randomize(population.begin(), population.end(), rng_);
231	0	std::vector<Candidate> shuffledPop1 = population;
232	0	randomize(population.begin(), population.end(), rng_);
233	0	std::vector<Candidate> shuffledPop2 = population;
234	0	randomize(population.begin(), population.end(), rng_);
235	0	mirrorPopulation = shuffledPop1;
236
237	0	adaptSizeWeights();
238
239	0	for (Size popIter = 0; popIter < population.size(); popIter++) {
240	0	if (rng_.nextReal() < 0.1){
241	0	population[popIter].values = rotateArray(bestMemberEver_.values);
242	0	}else {
243	0	population[popIter].values = bestMemberEver_.values
244	0	+ currGenSizeWeights_[popIter]
245	0	* (shuffledPop1[popIter].values - shuffledPop2[popIter].values);
246	0	}
247	0	}
248	0	}
249	0	break;
250
251	0	default:
252	0	QL_FAIL("Unknown strategy ("
253	0	<< Integer(configuration().strategy) << ")");
254	0	}
255		// in order to avoid unnecessary copying we use the same population object for mutants
256	0	crossover(oldPopulation, population, population, mirrorPopulation, p);
257	0	}
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	0	Problem& p) const {
265
266	0	if (configuration().crossoverIsAdaptive) {
267	0	adaptCrossover();
268	0	}
269
270	0	Array mutationProbabilities = getMutationProbabilities(population);
271
272	0	std::vector<Array> crossoverMask(population.size(),
273	0	Array(population.front().values.size(), 1.0));
274	0	std::vector<Array> invCrossoverMask = crossoverMask;
275	0	getCrossoverMask(crossoverMask, invCrossoverMask, mutationProbabilities);
276
277		// crossover of the old and mutant population
278	0	for (Size popIter = 0; popIter < population.size(); popIter++) {
279	0	population[popIter].values = oldPopulation[popIter].values * invCrossoverMask[popIter]
280	0	+ mutantPopulation[popIter].values * crossoverMask[popIter];
281		// immediately apply bounds if specified
282	0	if (configuration().applyBounds) {
283	0	for (Size memIter = 0; memIter < population[popIter].values.size(); memIter++) {
284	0	if (population[popIter].values[memIter] > upperBound_[memIter])
285	0	population[popIter].values[memIter] = upperBound_[memIter]
286	0	+ rng_.nextReal()
287	0	* (mirrorPopulation[popIter].values[memIter]
288	0	- upperBound_[memIter]);
289	0	if (population[popIter].values[memIter] < lowerBound_[memIter])
290	0	population[popIter].values[memIter] = lowerBound_[memIter]
291	0	+ rng_.nextReal()
292	0	* (mirrorPopulation[popIter].values[memIter]
293	0	- lowerBound_[memIter]);
294	0	}
295	0	}
296		// evaluate objective function as soon as possible to avoid unnecessary loops
297	0	try {
298	0	population[popIter].cost = p.value(population[popIter].values);
299	0	} catch (Error&) {
300	0	population[popIter].cost = QL_MAX_REAL;
301	0	}
302	0	if (!std::isfinite(population[popIter].cost))
303	0	population[popIter].cost = QL_MAX_REAL;
304
305	0	}
306	0	}
307
308		void DifferentialEvolution::getCrossoverMask(
309		std::vector<Array> & crossoverMask,
310		std::vector<Array> & invCrossoverMask,
311	0	const Array & mutationProbabilities) const {
312	0	for (Size cmIter = 0; cmIter < crossoverMask.size(); cmIter++) {
313	0	for (Size memIter = 0; memIter < crossoverMask[cmIter].size(); memIter++) {
314	0	if (rng_.nextReal() < mutationProbabilities[cmIter]) {
315	0	invCrossoverMask[cmIter][memIter] = 0.0;
316	0	} else {
317	0	crossoverMask[cmIter][memIter] = 0.0;
318	0	}
319	0	}
320	0	}
321	0	}
322
323		Array DifferentialEvolution::getMutationProbabilities(
324	0	const std::vector<Candidate> & population) const {
325	0	Array mutationProbabilities = currGenCrossover_;
326	0	switch (configuration().crossoverType) {
327	0	case Normal:
328	0	break;
329	0	case Binomial:
330	0	mutationProbabilities = currGenCrossover_
331	0	* (1.0 - 1.0 / population.front().values.size())
332	0	+ 1.0 / population.front().values.size();
333	0	break;
334	0	case Exponential:
335	0	for (Size coIter = 0;coIter< currGenCrossover_.size(); coIter++){
336	0	mutationProbabilities[coIter] =
337	0	(1.0 - std::pow(currGenCrossover_[coIter],
338	0	(int) population.front().values.size()))
339	0	/ (population.front().values.size()
340	0	* (1.0 - currGenCrossover_[coIter]));
341	0	}
342	0	break;
343	0	default:
344	0	QL_FAIL("Unknown crossover type ("
345	0	<< Integer(configuration().crossoverType) << ")");
346	0	break;
347	0	}
348	0	return mutationProbabilities;
349	0	}
350
351	0	Array DifferentialEvolution::rotateArray(Array a) const {
352	0	randomize(a.begin(), a.end(), rng_);
353	0	return a;
354	0	}
355
356	0	void DifferentialEvolution::adaptSizeWeights() const {
357		// [=Fl & =Fu] respectively see Brest, J. et al., 2006,
358		// "Self-Adapting Control Parameters in Differential
359		// Evolution"
360	0	Real sizeWeightLowerBound = 0.1, sizeWeightUpperBound = 0.9;
361		// [=tau1] A Comparative Study on Numerical Benchmark
362		// Problems." page 649 for reference
363	0	Real sizeWeightChangeProb = 0.1;
364	0	for (Real& currGenSizeWeight : currGenSizeWeights_) {
365	0	if (rng_.nextReal() < sizeWeightChangeProb)
366	0	currGenSizeWeight = sizeWeightLowerBound + rng_.nextReal() * sizeWeightUpperBound;
367	0	}
368	0	}
369
370	0	void DifferentialEvolution::adaptCrossover() const {
371	0	Real crossoverChangeProb = 0.1; // [=tau2]
372	0	for (Real& coIter : currGenCrossover_) {
373	0	if (rng_.nextReal() < crossoverChangeProb)
374	0	coIter = rng_.nextReal();
375	0	}
376	0	}
377
378		void DifferentialEvolution::fillInitialPopulation(
379		std::vector<Candidate> & population,
380	0	const Problem& p) const {
381
382		// use initial values provided by the user
383	0	population.front().values = p.currentValue();
384	0	population.front().cost = p.costFunction().value(population.front().values);
385		// rest of the initial population is random
386	0	for (Size j = 1; j < population.size(); ++j) {
387	0	for (Size i = 0; i < p.currentValue().size(); ++i) {
388	0	Real l = lowerBound_[i], u = upperBound_[i];
389	0	population[j].values[i] = l + (u-l)*rng_.nextReal();
390	0	}
391	0	population[j].cost = p.costFunction().value(population[j].values);
392	0	if (!std::isfinite(population[j].cost))
393	0	population[j].cost = QL_MAX_REAL;
394	0	}
395	0	}
396
397		}