QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
fireflyalgorithm.cpp
Go to the documentation of this file.
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4Copyright (C) 2015 Andres Hernandez
5
6This file is part of QuantLib, a free-software/open-source library
7for financial quantitative analysts and developers - http://quantlib.org/
8
9QuantLib is free software: you can redistribute it and/or modify it
10under the terms of the QuantLib license. You should have received a
11copy of the license along with this program; if not, please email
12<quantlib-dev@lists.sf.net>. The license is also available online at
13<http://quantlib.org/license.shtml>.
14
15This program is distributed in the hope that it will be useful, but WITHOUT
16ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
22#include <algorithm>
23#include <cmath>
24#include <utility>
25
26namespace QuantLib {
28 ext::shared_ptr<Intensity> intensity,
29 ext::shared_ptr<RandomWalk> randomWalk,
30 Size Mde,
31 Real mutation,
32 Real crossover,
33 unsigned long seed)
34 : mutation_(mutation), crossover_(crossover), M_(M), Mde_(Mde), Mfa_(M_ - Mde_),
35 intensity_(std::move(intensity)), randomWalk_(std::move(randomWalk)),
36 generator_(seed), distribution_(Mfa_, Mde > 0 ? M_ - 1 : M_),
37 rng_(seed) {
39 "Differential Evolution subpopulation cannot be larger than total population");
40 }
41
42 void FireflyAlgorithm::startState(Problem &P, const EndCriteria &endCriteria) {
43 N_ = P.currentValue().size();
44 x_.reserve(M_);
45 xI_.reserve(M_);
46 xRW_.reserve(M_);
47 values_.reserve(M_);
50 Array bounds = uX_ - lX_;
51
52 //Random initialization is done by Sobol sequence
53 SobolRsg sobol(N_);
54
55 //Prepare containers
56 for (Size i = 0; i < M_; i++) {
58 x_.emplace_back(N_, 0.0);
59 xI_.emplace_back(N_, 0.0);
60 xRW_.emplace_back(N_, 0.0);
61 Array& x = x_.back();
62 for (Size j = 0; j < N_; j++) {
63 //Assign X=lb+(ub-lb)*random
64 x[j] = lX_[j] + bounds[j] * sample[j];
65 }
66 //Evaluate point
67 values_.emplace_back(P.value(x), i);
68 }
69
70 //init intensity & randomWalk
71 intensity_->init(this);
72 randomWalk_->init(this);
73 }
74
76 QL_REQUIRE(!P.constraint().empty(), "Firefly Algorithm is a constrained optimizer");
78 P.reset();
79 Size iteration = 0;
80 Size iterationStat = 0;
81 Size maxIteration = endCriteria.maxIterations();
82 Size maxIStationary = endCriteria.maxStationaryStateIterations();
83
84 startState(P, endCriteria);
85
86 bool isFA = Mfa_ > 0;
87 //Variables for DE
88 Array z(N_, 0.0);
89 Size indexR1, indexR2;
90 decltype(distribution_)::param_type nParam(0, N_ - 1);
91
92 //Set best value & position
93 Real bestValue = values_[0].first;
94 Size bestPosition = 0;
95 for (Size i = 1; i < M_; i++) {
96 if (values_[i].first < bestValue) {
97 bestPosition = i;
98 bestValue = values_[i].first;
99 }
100 }
101 Array bestX = x_[bestPosition];
102
103 //Run optimization
104 do {
105 iteration++;
106 iterationStat++;
107 //Check if stopping criteria is met
108 if (iteration > maxIteration || iterationStat > maxIStationary)
109 break;
110
111 //Divide into two subpopulations
112 //First sort values
113 std::sort(values_.begin(), values_.end());
114
115 //Differential evolution
116 if(Mfa_ < M_){
117 Size indexBest = values_[0].second;
118 Array& xBest = x_[indexBest];
119 for (Size i = Mfa_; i < M_; i++) {
120 if (!isFA) {
121 //Pure DE requires random index
122 indexBest = distribution_(generator_);
123 xBest = x_[indexBest];
124 }
125 do {
126 indexR1 = distribution_(generator_);
127 } while(indexR1 == indexBest);
128 do {
129 indexR2 = distribution_(generator_);
130 } while(indexR2 == indexBest || indexR2 == indexR1);
131
132 Size index = values_[i].second;
133 Array& x = x_[index];
134 Array& xR1 = x_[indexR1];
135 Array& xR2 = x_[indexR2];
136 Size rIndex = distribution_(generator_, nParam);
137 for (Size j = 0; j < N_; j++) {
138 if (j == rIndex || rng_.nextReal() <= crossover_) {
139 //Change x[j] according to crossover
140 z[j] = xBest[j] + mutation_*(xR1[j] - xR2[j]);
141 } else {
142 z[j] = x[j];
143 }
144 //Enforce bounds on positions
145 if (z[j] < lX_[j]) {
146 z[j] = lX_[j];
147 }
148 else if (z[j] > uX_[j]) {
149 z[j] = uX_[j];
150 }
151 }
152 Real val = P.value(z);
153 if (val < values_[index].first) {
154 //Accept new point
155 x = z;
156 values_[index].first = val;
157 //mark best
158 if (val < bestValue) {
159 bestValue = val;
160 bestX = x;
161 iterationStat = 0;
162 }
163 }
164 }
165 }
166
167 //Firefly algorithm
168 if(isFA){
169 //According to the intensity, determine best global position
170 intensity_->findBrightest();
171
172 //Prepare random walk
173 randomWalk_->walk();
174
175 //Loop over particles
176 for (Size i = 0; i < Mfa_; i++) {
177 Size index = values_[i].second;
178 Array& x = x_[index];
179 Array& xI = xI_[index];
180 Array& xRW = xRW_[index];
181
182 //Loop over dimensions
183 for (Size j = 0; j < N_; j++) {
184 //Update position
185 z[j] = x[j] + xI[j] + xRW[j];
186 //Enforce bounds on positions
187 if (z[j] < lX_[j]) {
188 z[j] = lX_[j];
189 }
190 else if (z[j] > uX_[j]) {
191 z[j] = uX_[j];
192 }
193 }
194 Real val = P.value(z);
195 if(!std::isnan(val))
196 {
197 //Accept new point
198 x = z;
199 values_[index].first = val;
200 //mark best
201 if (val < bestValue) {
202 bestValue = val;
203 bestX = x;
204 iterationStat = 0;
205 }
206 }
207 }
208 }
209 } while (true);
210 if (iteration > maxIteration)
212 else
214
215 //Set result to best point
216 P.setCurrentValue(bestX);
217 P.setFunctionValue(bestValue);
218 return ecType;
219 }
220
222 //Brightest ignores all others
223 Array& xI = (*xI_)[(*values_)[0].second];
224 for (Size j = 0; j < N_; j++) {
225 xI[j] = 0.0;
226 }
227
228 for (Size i = 1; i < Mfa_; i++) {
229 //values_ is already sorted
230 Size index = (*values_)[i].second;
231 const Array& x = (*x_)[index];
232 Array& xI = (*xI_)[index];
233 for (Size j = 0; j < N_; j++) {
234 xI[j] = 0.0;
235 }
236 Real valueX = (*values_)[i].first;
237 for (Size k = 0; k < i - 1; k++){
238 const Array& y = (*x_)[(*values_)[k].second];
239 Real valueY = (*values_)[k].first;
240 Real intensity = intensityImpl(valueX, valueY, distance(x, y));
241 xI += intensity*(y - x);
242 }
243 }
244 }
245}
246
1-D array used in linear algebra.
Definition: array.hpp:52
Size size() const
dimension of the array
Definition: array.hpp:495
bool empty() const
Definition: constraint.hpp:56
Array lowerBound(const Array &params) const
Definition: constraint.hpp:66
Array upperBound(const Array &params) const
Definition: constraint.hpp:58
Criteria to end optimization process:
Definition: endcriteria.hpp:40
Size maxIterations() const
Size maxStationaryStateIterations() const
virtual Real intensityImpl(Real valueX, Real valueY, Real distance)=0
Real distance(const Array &x, const Array &y) const
void findBrightest()
find brightest firefly for each firefly
std::vector< std::pair< Real, Size > > values_
std::uniform_int_distribution< QuantLib::Size > distribution_
ext::shared_ptr< Intensity > intensity_
void startState(Problem &P, const EndCriteria &endCriteria)
MersenneTwisterUniformRng rng_
std::vector< Array > xRW_
FireflyAlgorithm(Size M, ext::shared_ptr< Intensity > intensity, ext::shared_ptr< RandomWalk > randomWalk, Size Mde=0, Real mutationFactor=1.0, Real crossoverFactor=0.5, unsigned long seed=SeedGenerator::instance().get())
std::vector< Array > xI_
EndCriteria::Type minimize(Problem &P, const EndCriteria &endCriteria) override
minimize the optimization problem P
ext::shared_ptr< RandomWalk > randomWalk_
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
void setCurrentValue(const Array &currentValue)
Definition: problem.hpp:76
Sobol low-discrepancy sequence generator.
Definition: sobolrsg.hpp:110
const SobolRsg::sample_type & nextSequence() const
Definition: sobolrsg.hpp:134
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
Implementation based on: Yang, Xin-She (2009) Firefly Algorithm, Levy Flights and Global Optimization...
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
Definition: any.hpp:35
STL namespace.
Sobol low-discrepancy sequence generator.