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/particleswarmoptimization.hpp>
21#include <ql/math/randomnumbers/sobolrsg.hpp>
22#include <cmath>
23#include <utility>
24
25using std::sqrt;
26
27namespace QuantLib {
28 ParticleSwarmOptimization::ParticleSwarmOptimization(Size M,
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(x: 2.0 - phi - sqrt(x: phi*phi - 4 * phi));
38 c1_ = c0_*c1;
39 c2_ = c0_*c2;
40 }
41
42 ParticleSwarmOptimization::ParticleSwarmOptimization(Size M,
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
52 void ParticleSwarmOptimization::startState(Problem &P, const EndCriteria &endCriteria) {
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: M_, N: N_, c0: c0_, endCriteria);
58 X_.reserve(n: M_);
59 V_.reserve(n: M_);
60 pBX_.reserve(n: M_);
61 pBF_ = Array(M_);
62 gBX_.reserve(n: M_);
63 gBF_ = Array(M_);
64 uX_ = P.constraint().upperBound(params: P.currentValue());
65 lX_ = P.constraint().lowerBound(params: P.currentValue());
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++) {
73 const SobolRsg::sample_type::value_type &sample = sobol.nextSequence().value;
74 X_.emplace_back(args&: N_, args: 0.0);
75 Array& x = X_.back();
76 V_.emplace_back(args&: N_, args: 0.0);
77 Array& v = V_.back();
78 gBX_.emplace_back(args&: N_, args: 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: X_.back());
87 pBF_[i] = P.value(x: X_.back());
88 }
89
90 //init topology & inertia
91 topology_->init(pso: this);
92 inertia_->init(pso: this);
93 }
94
95 EndCriteria::Type ParticleSwarmOptimization::minimize(Problem &P, const EndCriteria &endCriteria) {
96 QL_REQUIRE(!P.constraint().empty(), "PSO is a constrained optimizer");
97
98 EndCriteria::Type ecType = EndCriteria::None;
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)
169 ecType = EndCriteria::MaxIterations;
170 else
171 ecType = EndCriteria::StationaryPoint;
172
173 //Set result to best point
174 P.setCurrentValue(pBX_[bestPosition]);
175 P.setFunctionValue(bestValue);
176 return ecType;
177 }
178
179 void AdaptiveInertia::setValues() {
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;
187 adaptiveCounter--;
188 }
189 else {
190 adaptiveCounter++;
191 }
192 if (adaptiveCounter > sh_) {
193 c0_ = std::max(a: minInertia_, b: std::min(a: maxInertia_, b: c0_*0.5));
194 }
195 else if (adaptiveCounter < sl_) {
196 c0_ = std::max(a: minInertia_, b: std::min(a: maxInertia_, b: 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
208 void KNeighbors::findSocialBest() {
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(a: i + K_, b: M_);
214 //Search K_ neighbors downwards
215 Size lower = std::max(a: i, b: 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
243 ClubsTopology::ClubsTopology(Size defaultClubs,
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_) {
252 QL_REQUIRE(totalClubs_ >= defaultClubs_,
253 "Total number of clubs must be larger or equal than default clubs");
254 QL_REQUIRE(defaultClubs_ >= minClubs_,
255 "Number of default clubs must be larger or equal than minimum clubs");
256 QL_REQUIRE(maxClubs_ >= defaultClubs_,
257 "Number of maximum clubs must be larger or equal than default clubs");
258 QL_REQUIRE(totalClubs_ >= maxClubs_,
259 "Total number of clubs must be larger or equal than maximum clubs");
260 }
261
262 void ClubsTopology::setSize(Size M) {
263 M_ = M;
264
265 if (defaultClubs_ < totalClubs_) {
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++) {
272 Size index = distribution_(generator_);
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
286 void ClubsTopology::findSocialBest() {
287 //Update iteration
288 iteration_++;
289 bool reset = false;
290 if (iteration_ == resetIteration_) {
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(particle: i, currentClubs);
338 }
339 else if (worst) {
340 //Join random club
341 joinRandomClub(particle: 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(particle: i, currentClubs);
348 }
349 else {
350 //Leave random club
351 leaveRandomClub(particle: 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
403

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