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#include <ql/experimental/math/fireflyalgorithm.hpp>
21#include <ql/math/randomnumbers/sobolrsg.hpp>
22#include <algorithm>
23#include <cmath>
24#include <utility>
25
26namespace QuantLib {
27 FireflyAlgorithm::FireflyAlgorithm(Size M,
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) {
38 QL_REQUIRE(M_ >= Mde_,
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(n: M_);
45 xI_.reserve(n: M_);
46 xRW_.reserve(n: M_);
47 values_.reserve(n: M_);
48 uX_ = P.constraint().upperBound(params: P.currentValue());
49 lX_ = P.constraint().lowerBound(params: P.currentValue());
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++) {
57 const SobolRsg::sample_type::value_type &sample = sobol.nextSequence().value;
58 x_.emplace_back(args&: N_, args: 0.0);
59 xI_.emplace_back(args&: N_, args: 0.0);
60 xRW_.emplace_back(args&: N_, args: 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(args: P.value(x), args&: i);
68 }
69
70 //init intensity & randomWalk
71 intensity_->init(fa: this);
72 randomWalk_->init(fa: this);
73 }
74
75 EndCriteria::Type FireflyAlgorithm::minimize(Problem &P, const EndCriteria &endCriteria) {
76 QL_REQUIRE(!P.constraint().empty(), "Firefly Algorithm is a constrained optimizer");
77 EndCriteria::Type ecType = EndCriteria::None;
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(first: values_.begin(), last: 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(x: 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(x: z);
195 if(!std::isnan(x: 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)
211 ecType = EndCriteria::MaxIterations;
212 else
213 ecType = EndCriteria::StationaryPoint;
214
215 //Set result to best point
216 P.setCurrentValue(bestX);
217 P.setFunctionValue(bestValue);
218 return ecType;
219 }
220
221 void FireflyAlgorithm::Intensity::findBrightest() {
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: distance(x, y));
241 xI += intensity*(y - x);
242 }
243 }
244 }
245}
246
247

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