QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
particleswarmoptimization.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 <cmath>
23#include <utility>
24
25using std::sqrt;
26
27namespace QuantLib {
29 ext::shared_ptr<Topology> topology,
30 ext::shared_ptr<Inertia> inertia,
31 Real c1,
32 Real c2,
33 unsigned long seed)
34 : M_(M), rng_(seed), topology_(std::move(topology)), inertia_(std::move(inertia)) {
35 Real phi = c1 + c2;
36 QL_ENSURE(phi*phi - 4 * phi != 0.0, "Invalid phi");
37 c0_ = 2.0 / std::abs(2.0 - phi - sqrt(phi*phi - 4 * phi));
38 c1_ = c0_*c1;
39 c2_ = c0_*c2;
40 }
41
43 ext::shared_ptr<Topology> topology,
44 ext::shared_ptr<Inertia> inertia,
45 Real omega,
46 Real c1,
47 Real c2,
48 unsigned long seed)
49 : M_(M), c0_(omega), c1_(c1), c2_(c2), rng_(seed), topology_(std::move(topology)),
50 inertia_(std::move(inertia)) {}
51
53 QL_REQUIRE(topology_, "Invalid topology");
54 QL_REQUIRE(inertia_, "Invalid inertia");
55 N_ = P.currentValue().size();
56 topology_->setSize(M_);
57 inertia_->setSize(M_, N_, c0_, endCriteria);
58 X_.reserve(M_);
59 V_.reserve(M_);
60 pBX_.reserve(M_);
61 pBF_ = Array(M_);
62 gBX_.reserve(M_);
63 gBF_ = Array(M_);
66 Array bounds = uX_ - lX_;
67
68 //Random initialization is done by Sobol sequence
69 SobolRsg sobol(N_ * 2);
70
71 //Prepare containers
72 for (Size i = 0; i < M_; i++) {
74 X_.emplace_back(N_, 0.0);
75 Array& x = X_.back();
76 V_.emplace_back(N_, 0.0);
77 Array& v = V_.back();
78 gBX_.emplace_back(N_, 0.0);
79 for (Size j = 0; j < N_; j++) {
80 //Assign X=lb+(ub-lb)*random
81 x[j] = lX_[j] + bounds[j] * sample[2 * j];
82 //Assign V=(ub-lb)*2*random-(ub-lb) -> between (lb-ub) and (ub-lb)
83 v[j] = bounds[j] * (2.0*sample[2 * j + 1] - 1.0);
84 }
85 //Evaluate X and assign as personal best
86 pBX_.push_back(X_.back());
87 pBF_[i] = P.value(X_.back());
88 }
89
90 //init topology & inertia
91 topology_->init(this);
92 inertia_->init(this);
93 }
94
96 QL_REQUIRE(!P.constraint().empty(), "PSO is a constrained optimizer");
97
99 P.reset();
100 Size iteration = 0;
101 Size iterationStat = 0;
102 Size maxIteration = endCriteria.maxIterations();
103 Size maxIStationary = endCriteria.maxStationaryStateIterations();
104 Real bestValue = QL_MAX_REAL;
105 Size bestPosition = 0;
106
107 startState(P, endCriteria);
108 //Set best value & position
109 for (Size i = 0; i < M_; i++) {
110 if (pBF_[i] < bestValue) {
111 bestValue = pBF_[i];
112 bestPosition = i;
113 }
114 }
115
116 //Run optimization
117 do {
118 iteration++;
119 iterationStat++;
120 //Check if stopping criteria is met
121 if (iteration > maxIteration || iterationStat > maxIStationary)
122 break;
123
124 //According to the topology, determine best global position
125 topology_->findSocialBest();
126
127 //Call inertia to change internal state
128 inertia_->setValues();
129
130 //Loop over particles
131 for (Size i = 0; i < M_; i++) {
132 Array& x = X_[i];
133 Array& pB = pBX_[i];
134 const Array& gB = gBX_[i];
135 Array& v = V_[i];
136
137 //Loop over dimensions
138 for (Size j = 0; j < N_; j++) {
139 //Update velocity
140 v[j] += c1_*rng_.nextReal()*(pB[j] - x[j]) + c2_*rng_.nextReal()*(gB[j] - x[j]);
141 //Update position
142 x[j] += v[j];
143 //Enforce bounds on positions
144 if (x[j] < lX_[j]) {
145 x[j] = lX_[j];
146 v[j] = 0.0;
147 }
148 else if (x[j] > uX_[j]) {
149 x[j] = uX_[j];
150 v[j] = 0.0;
151 }
152 }
153 //Evaluate x
154 Real f = P.value(x);
155 if (f < pBF_[i]) {
156 //Update personal best
157 pBF_[i] = f;
158 pB = x;
159 //Check stationary condition
160 if (f < bestValue) {
161 bestValue = f;
162 bestPosition = i;
163 iterationStat = 0;
164 }
165 }
166 }
167 } while (true);
168 if (iteration > maxIteration)
170 else
172
173 //Set result to best point
174 P.setCurrentValue(pBX_[bestPosition]);
175 P.setFunctionValue(bestValue);
176 return ecType;
177 }
178
180 Real currBest = (*pBF_)[0];
181 for (Size i = 1; i < M_; i++) {
182 if (currBest >(*pBF_)[i]) currBest = (*pBF_)[i];
183 }
184 if (started_) { //First iteration leaves inertia unchanged
185 if (currBest < best_) {
186 best_ = currBest;
188 }
189 else {
191 }
192 if (adaptiveCounter > sh_) {
193 c0_ = std::max(minInertia_, std::min(maxInertia_, c0_*0.5));
194 }
195 else if (adaptiveCounter < sl_) {
196 c0_ = std::max(minInertia_, std::min(maxInertia_, c0_*2.0));
197 }
198 }
199 else {
200 best_ = currBest;
201 started_ = true;
202 }
203 for (Size i = 0; i < M_; i++) {
204 (*V_)[i] *= c0_;
205 }
206 }
207
209 for (Size i = 0; i < M_; i++) {
210 Real bestF = (*pBF_)[i];
211 Size bestX = 0;
212 //Search K_ neightbors upwards
213 Size upper = std::min(i + K_, M_);
214 //Search K_ neighbors downwards
215 Size lower = std::max(i, K_ + 1) - K_ - 1;
216 for (Size j = lower; j < upper; j++) {
217 if ((*pBF_)[j] < bestF) {
218 bestF = (*pBF_)[j];
219 bestX = j;
220 }
221 }
222 if (i + K_ >= M_) { //loop around if i+K >= M_
223 for (Size j = 0; j < i + K_ - M_; j++) {
224 if ((*pBF_)[j] < bestF) {
225 bestF = (*pBF_)[j];
226 bestX = j;
227 }
228 }
229 }
230 else if (i < K_) {//loop around from above
231 for (Size j = M_ - (K_ - i) - 1; j < M_; j++) {
232 if ((*pBF_)[j] < bestF) {
233 bestF = (*pBF_)[j];
234 bestX = j;
235 }
236 }
237 }
238 (*gBX_)[i] = (*pBX_)[bestX];
239 (*gBF_)[i] = bestF;
240 }
241 }
242
244 Size totalClubs,
245 Size maxClubs,
246 Size minClubs,
247 Size resetIteration,
248 unsigned long seed)
249 : totalClubs_(totalClubs), maxClubs_(maxClubs), minClubs_(minClubs),
250 defaultClubs_(defaultClubs), resetIteration_(resetIteration), bestByClub_(totalClubs, 0),
251 worstByClub_(totalClubs, 0), generator_(seed), distribution_(1, totalClubs_) {
253 "Total number of clubs must be larger or equal than default clubs");
255 "Number of default clubs must be larger or equal than minimum clubs");
257 "Number of maximum clubs must be larger or equal than default clubs");
259 "Total number of clubs must be larger or equal than maximum clubs");
260 }
261
263 M_ = M;
264
266 clubs4particles_ = std::vector<std::vector<bool> >(M_, std::vector<bool>(totalClubs_, false));
267 particles4clubs_ = std::vector<std::vector<bool> >(totalClubs_, std::vector<bool>(M_, false));
268 //Assign particles to clubs randomly
269 for (Size i = 0; i < M_; i++) {
270 std::vector<bool> &clubSet = clubs4particles_[i];
271 for (Size j = 0; j < defaultClubs_; j++) {
273 while (clubSet[index]) { index = distribution_(generator_); }
274 clubSet[index] = true;
275 particles4clubs_[index][i] = true;
276 }
277 }
278 }
279 else {
280 //Since totalClubs_ == defaultClubs_, then just initialize to true
281 clubs4particles_ = std::vector<std::vector<bool> >(M_, std::vector<bool>(totalClubs_, true));
282 particles4clubs_ = std::vector<std::vector<bool> >(totalClubs_, std::vector<bool>(M_, true));
283 }
284 }
285
287 //Update iteration
288 iteration_++;
289 bool reset = false;
291 iteration_ = 0;
292 reset = true;
293 }
294
295 //Find best by current club
296 for (Size i = 0; i < totalClubs_; i++) {
297 Real bestByClub = QL_MAX_REAL;
298 Real worstByClub = -QL_MAX_REAL;
299 Size bestP = 0;
300 Size worstP = 0;
301 const std::vector<bool> &particlesSet = particles4clubs_[i];
302 for (Size j = 0; j < M_; j++) {
303 if (particlesSet[j]) {
304 if (bestByClub >(*pBF_)[j]) {
305 bestByClub = (*pBF_)[j];
306 bestP = j;
307 }
308 else if (worstByClub < (*pBF_)[j]) {
309 worstByClub = (*pBF_)[j];
310 worstP = j;
311 }
312 }
313 }
314 bestByClub_[i] = bestP;
315 worstByClub_[i] = worstP;
316 }
317
318 //Update clubs && global best
319 for (Size i = 0; i < M_; i++) {
320 std::vector<bool> &clubSet = clubs4particles_[i];
321 bool best = true;
322 bool worst = true;
323 Size currentClubs = 0;
324 for (Size j = 0; j < totalClubs_; j++) {
325 if (clubSet[j]) {
326 //If still thought of the best, check if best in club j
327 if (best && i != bestByClub_[j]) best = false;
328 //If still thought of the worst, check if worst in club j
329 if (worst && i != worstByClub_[j]) worst = false;
330 //Update currentClubs
331 currentClubs++;
332 }
333 }
334 //Update clubs
335 if (best) {
336 //Leave random club
337 leaveRandomClub(i, currentClubs);
338 }
339 else if (worst) {
340 //Join random club
341 joinRandomClub(i, currentClubs);
342 }
343 else if (reset && currentClubs != defaultClubs_) {
344 //If membership != defaultClubs_, then leave or join accordingly
345 if (currentClubs < defaultClubs_) {
346 //Join random club
347 joinRandomClub(i, currentClubs);
348 }
349 else {
350 //Leave random club
351 leaveRandomClub(i, currentClubs);
352 }
353 }
354
355 //Update global best
356 Real bestNeighborF = QL_MAX_REAL;
357 Size bestNeighborX = 0;
358 for (Size j = 0; j < totalClubs_; j++) {
359 if (clubSet[j] && bestNeighborF >(*pBF_)[bestByClub_[j]]) {
360 bestNeighborF = (*pBF_)[bestByClub_[j]];
361 bestNeighborX = j;
362 }
363 }
364 (*gBX_)[i] = (*pBX_)[bestNeighborX];
365 (*gBF_)[i] = bestNeighborF;
366 }
367 }
368
369 void ClubsTopology::leaveRandomClub(Size particle, Size currentClubs) {
370 Size randIndex = distribution_(generator_, param_type(1, currentClubs));
371 Size index = 1;
372 std::vector<bool> &clubSet = clubs4particles_[particle];
373 for (Size j = 0; j < totalClubs_; j++) {
374 if (clubSet[j]) {
375 if (index == randIndex) {
376 clubSet[j] = false;
377 particles4clubs_[j][particle] = false;
378 break;
379 }
380 index++;
381 }
382 }
383 }
384
385 void ClubsTopology::joinRandomClub(Size particle, Size currentClubs) {
386 Size randIndex = totalClubs_ == currentClubs ? 1 :
387 distribution_(generator_, param_type(1, totalClubs_ - currentClubs));
388 Size index = 1;
389 std::vector<bool> &clubSet = clubs4particles_[particle];
390 for (Size j = 0; j < totalClubs_; j++) {
391 if (!clubSet[j]) {
392 if (index == randIndex) {
393 clubSet[j] = true;
394 particles4clubs_[j][particle] = true;
395 break;
396 }
397 index++;
398 }
399 }
400 }
401}
402
void setValues() override
produce changes to PSO state for current iteration
1-D array used in linear algebra.
Definition: array.hpp:52
Size size() const
dimension of the array
Definition: array.hpp:495
std::uniform_int_distribution< QuantLib::Size > distribution_
std::vector< std::vector< bool > > particles4clubs_
std::vector< std::vector< bool > > clubs4particles_
ClubsTopology(Size defaultClubs, Size totalClubs, Size maxClubs, Size minClubs, Size resetIteration, unsigned long seed=SeedGenerator::instance().get())
void setSize(Size M) override
initialize state for current problem
decltype(distribution_)::param_type param_type
void findSocialBest() override
produce changes to PSO state for current iteration
void joinRandomClub(Size particle, Size currentClubs)
void leaveRandomClub(Size particle, Size currentClubs)
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
void findSocialBest() override
produce changes to PSO state for current iteration
Real nextReal() const
return a random number in the (0.0, 1.0)-interval
ParticleSwarmOptimization(Size M, ext::shared_ptr< Topology > topology, ext::shared_ptr< Inertia > inertia, Real c1=2.05, Real c2=2.05, unsigned long seed=SeedGenerator::instance().get())
void startState(Problem &P, const EndCriteria &endCriteria)
EndCriteria::Type minimize(Problem &P, const EndCriteria &endCriteria) override
minimize the optimization problem P
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_ENSURE(condition, message)
throw an error if the given post-condition is not verified
Definition: errors.hpp:130
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
#define QL_MAX_REAL
Definition: qldefines.hpp:176
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.
Implementation based on: Clerc, M., Kennedy, J. (2002) The particle swarm-explosion,...
ext::shared_ptr< BlackVolTermStructure > v
Sobol low-discrepancy sequence generator.