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
20/*! \file fireflyalgorithm.hpp
21\brief Implementation based on:
22Yang, Xin-She (2009) Firefly Algorithm, Levy Flights and Global
23Optimization. Research and Development in Intelligent Systems XXVI, pp 209-218.
24http://arxiv.org/pdf/1003.1464.pdf
25*/
26
27#ifndef quantlib_optimization_fireflyalgorithm_hpp
28#define quantlib_optimization_fireflyalgorithm_hpp
29
30#include <ql/math/optimization/problem.hpp>
31#include <ql/math/optimization/constraint.hpp>
32#include <ql/experimental/math/isotropicrandomwalk.hpp>
33#include <ql/experimental/math/levyflightdistribution.hpp>
34#include <ql/math/randomnumbers/mt19937uniformrng.hpp>
35#include <ql/math/randomnumbers/seedgenerator.hpp>
36
37#include <cmath>
38#include <random>
39
40namespace QuantLib {
41
42 /*! The main process is as follows:
43 M individuals are used to explore the N-dimensional parameter space:
44 \f$ X_{i}^k = (X_{i, 1}^k, X_{i, 2}^k, \ldots, X_{i, N}^k) \f$ is the kth-iteration
45 for the ith-individual. X is updated via the rule
46 \f[
47 X_{i, j}^{k+1} = X_{i, j}^k + I(X^k)_{i,j} + RandomWalk_{i,j}^k
48 \f]
49
50 The intensity function I(X) should be monotonic
51 The optimization stops either because the number of iterations has been reached
52 or because the stationary function value limit has been reached.
53
54 The current implementation extends the normal Firefly Algorithm with a
55 differential evolution (DE) optimizer according to:
56 Afnizanfaizal Abdullah, et al. "A New Hybrid Firefly Algorithm for Complex and
57 Nonlinear Problem". Volume 151 of the series Advances in Intelligent and Soft
58 Computing pp 673-680, 2012.
59 http://link.springer.com/chapter/10.1007%2F978-3-642-28765-7_81
60
61 In effect this implementation provides a fully fledged DE global optimizer
62 as well. The Firefly Algorithm was easy to combine with DE because it already
63 contained a step where the current solutions are sorted. The population is
64 then divided into two subpopulations based on their order. The subpopulation
65 with the best results are updated via the firefly algorithm. The worse
66 subpopulation is updated via the DE operator:
67 \f[
68 Y^{k+1} = X_{best}^k + F(X_{r1}^k - X_{r2}^k)
69 \f]
70 and
71 \f[
72 X_{i,j}^{k+1} = Y_{i,j}^{k+1}\ \text{if} R_{i,j} <= C
73 \f]
74 \f[
75 X_{i,j}^{k+1} = X_{i,j}^{k+1}\ \text{otherwise}
76 \f]
77 where C is the crossover constant, and R is a random uniformly distributed
78 number.
79 */
80 class FireflyAlgorithm : public OptimizationMethod {
81 public:
82 class RandomWalk;
83 class Intensity;
84 FireflyAlgorithm(Size M,
85 ext::shared_ptr<Intensity> intensity,
86 ext::shared_ptr<RandomWalk> randomWalk,
87 Size Mde = 0,
88 Real mutationFactor = 1.0,
89 Real crossoverFactor = 0.5,
90 unsigned long seed = SeedGenerator::instance().get());
91 void startState(Problem &P, const EndCriteria &endCriteria);
92 EndCriteria::Type minimize(Problem& P, const EndCriteria& endCriteria) override;
93
94 protected:
95 std::vector<Array> x_, xI_, xRW_;
96 std::vector<std::pair<Real, Size> > values_;
97 Array lX_, uX_;
98 Real mutation_, crossover_;
99 Size M_, N_, Mde_, Mfa_;
100 ext::shared_ptr<Intensity> intensity_;
101 ext::shared_ptr<RandomWalk> randomWalk_;
102 std::mt19937 generator_;
103 std::uniform_int_distribution<QuantLib::Size> distribution_;
104 MersenneTwisterUniformRng rng_;
105 };
106
107 //! Base intensity class
108 /*! Derived classes need to implement only intensityImpl
109 */
110 class FireflyAlgorithm::Intensity {
111 friend class FireflyAlgorithm;
112 public:
113 virtual ~Intensity() = default;
114 //! find brightest firefly for each firefly
115 void findBrightest();
116 protected:
117 Size Mfa_, N_;
118 const std::vector<Array> *x_;
119 const std::vector<std::pair<Real, Size> > *values_;
120 std::vector<Array> *xI_;
121
122 virtual Real intensityImpl(Real valueX, Real valueY, Real distance) = 0;
123 inline Real distance(const Array& x, const Array& y) const {
124 Real d = 0.0;
125 for (Size i = 0; i < N_; i++) {
126 Real diff = x[i] - y[i];
127 d += diff*diff;
128 }
129 return d;
130 }
131
132 private:
133 void init(FireflyAlgorithm *fa) {
134 x_ = &fa->x_;
135 xI_ = &fa->xI_;
136 values_ = &fa->values_;
137 Mfa_ = fa->Mfa_;
138 N_ = fa->N_;
139 }
140 };
141
142 //! Exponential Intensity
143 /* Exponentially decreasing intensity
144 */
145 class ExponentialIntensity : public FireflyAlgorithm::Intensity {
146 public:
147 ExponentialIntensity(Real beta0, Real betaMin, Real gamma)
148 : beta0_(beta0), betaMin_(betaMin), gamma_(gamma) {}
149 protected:
150 Real intensityImpl(Real valueX, Real valueY, Real d) override {
151 return (beta0_ - betaMin_) * std::exp(x: -gamma_ * d) + betaMin_;
152 }
153 Real beta0_, betaMin_, gamma_;
154 };
155
156 //! Inverse Square Intensity
157 /* Inverse law square
158 */
159 class InverseLawSquareIntensity : public FireflyAlgorithm::Intensity {
160 public:
161 InverseLawSquareIntensity(Real beta0, Real betaMin)
162 : beta0_(beta0), betaMin_(betaMin) {}
163 protected:
164 Real intensityImpl(Real valueX, Real valueY, Real d) override {
165 return (beta0_ - betaMin_) / (d + QL_EPSILON) + betaMin_;
166 }
167 Real beta0_, betaMin_;
168 };
169
170 //! Base Random Walk class
171 /*! Derived classes need to implement only walkImpl
172 */
173 class FireflyAlgorithm::RandomWalk {
174 friend class FireflyAlgorithm;
175 public:
176 virtual ~RandomWalk() = default;
177 //! perform random walk
178 void walk() {
179 for (Size i = 0; i < Mfa_; i++) {
180 walkImpl(xRW&: (*xRW_)[(*values_)[i].second]);
181 }
182 }
183 protected:
184 Size Mfa_, N_;
185 const std::vector<Array> *x_;
186 const std::vector<std::pair<Real, Size> > *values_;
187 std::vector<Array> *xRW_;
188 Array *lX_, *uX_;
189
190 virtual void walkImpl(Array & xRW) = 0;
191 virtual void init(FireflyAlgorithm *fa) {
192 x_ = &fa->x_;
193 xRW_ = &fa->xRW_;
194 values_ = &fa->values_;
195 Mfa_ = fa->Mfa_;
196 N_ = fa->N_;
197 lX_ = &fa->lX_;
198 uX_ = &fa->uX_;
199 }
200 };
201
202 //! Distribution Walk
203 /* Random walk given by distribution template parameter. The
204 distribution must be compatible with std::mt19937.
205 */
206 template <class Distribution>
207 class DistributionRandomWalk : public FireflyAlgorithm::RandomWalk {
208 public:
209 explicit DistributionRandomWalk(Distribution dist,
210 Real delta = 0.9,
211 unsigned long seed = SeedGenerator::instance().get())
212 : walkRandom_(std::mt19937(seed), std::move(dist), 1, Array(1, 1.0), seed),
213 delta_(delta) {}
214 protected:
215 void walkImpl(Array& xRW) override {
216 walkRandom_.nextReal(&xRW[0]);
217 xRW *= delta_;
218 }
219 void init(FireflyAlgorithm* fa) override {
220 FireflyAlgorithm::RandomWalk::init(fa);
221 walkRandom_.setDimension(N_, *lX_, *uX_);
222 }
223 IsotropicRandomWalk<Distribution, std::mt19937> walkRandom_;
224 Real delta_;
225 };
226
227 //! Gaussian Walk
228 /* Gaussian random walk
229 */
230 class GaussianWalk : public DistributionRandomWalk<std::normal_distribution<QuantLib::Real>> {
231 public:
232 explicit GaussianWalk(Real sigma,
233 Real delta = 0.9,
234 unsigned long seed = SeedGenerator::instance().get())
235 : DistributionRandomWalk<std::normal_distribution<QuantLib::Real>>(
236 std::normal_distribution<QuantLib::Real>(0.0, sigma), delta, seed){}
237 };
238
239 //! Levy Flight Random Walk
240 /* Levy flight random walk
241 */
242 class LevyFlightWalk : public DistributionRandomWalk<LevyFlightDistribution> {
243 public:
244 explicit LevyFlightWalk(Real alpha,
245 Real xm = 0.5,
246 Real delta = 0.9,
247 unsigned long seed = SeedGenerator::instance().get())
248 : DistributionRandomWalk<LevyFlightDistribution>(
249 LevyFlightDistribution(xm, alpha), delta, seed) {}
250 };
251
252 //! Decreasing Random Walk
253 /* Gaussian random walk, but size of step decreases with each iteration step
254 */
255 class DecreasingGaussianWalk : public GaussianWalk {
256 public:
257 explicit DecreasingGaussianWalk(
258 Real sigma,
259 Real delta = 0.9,
260 unsigned long seed = SeedGenerator::instance().get())
261 : GaussianWalk(sigma, delta, seed), delta0_(delta) {}
262 protected:
263 void walkImpl(Array& xRW) override {
264 iteration_++;
265 if (iteration_ > Mfa_) {
266 //Every time all the fireflies have been processed
267 //multiply delta by itself
268 iteration_ = 0;
269 delta_ *= delta_;
270 }
271 GaussianWalk::walkImpl(xRW);
272 }
273 void init(FireflyAlgorithm* fa) override {
274 GaussianWalk::init(fa);
275 iteration_ = 0;
276 delta_ = delta0_;
277 }
278
279 private:
280 Real delta0_;
281 Size iteration_;
282 };
283}
284
285#endif
286

source code of quantlib/ql/experimental/math/fireflyalgorithm.hpp