1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*!
4 Copyright (C) 2016 Andres Hernandez
5
6 This file is part of QuantLib, a free-software/open-source library
7 for financial quantitative analysts and developers - http://quantlib.org/
8
9 QuantLib is free software: you can redistribute it and/or modify it
10 under the terms of the QuantLib license. You should have received a
11 copy 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
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
20#include <ql/qldefines.hpp>
21#if !defined(BOOST_ALL_NO_LIB) && defined(BOOST_MSVC)
22# include <ql/auto_link.hpp>
23#endif
24#include <ql/experimental/math/fireflyalgorithm.hpp>
25#include <ql/experimental/math/hybridsimulatedannealing.hpp>
26#include <ql/experimental/math/particleswarmoptimization.hpp>
27#include <ql/functional.hpp>
28#include <ql/math/optimization/differentialevolution.hpp>
29#include <ql/math/optimization/simulatedannealing.hpp>
30#include <ql/tuple.hpp>
31#include <iomanip>
32#include <iostream>
33#include <utility>
34
35using namespace QuantLib;
36
37unsigned long seed = 127;
38
39/*
40 Some benchmark functions taken from
41 https://en.wikipedia.org/wiki/Test_functions_for_optimization
42
43 Global optimizers have generally a lot of hyper-parameters, and one
44 * usually requires some hyper-parameter optimization to find appropriate values
45*/
46
47Real ackley(const Array& x) {
48 //Minimum is found at 0
49 Real p1 = 0.0, p2 = 0.0;
50
51 for (Real i : x) {
52 p1 += i * i;
53 p2 += std::cos(M_TWOPI * i);
54 }
55 p1 = -0.2*std::sqrt(x: 0.5*p1);
56 p2 *= 0.5;
57 return M_E + 20.0 - 20.0*std::exp(x: p1)-std::exp(x: p2);
58}
59
60Array ackleyValues(const Array& x) {
61 Array y(x.size());
62 for (Size i = 0; i < x.size(); i++) {
63 Real p1 = x[i] * x[i];
64 p1 = -0.2*std::sqrt(x: 0.5*p1);
65 Real p2 = 0.5*std::cos(M_TWOPI*x[i]);
66 y[i] = M_E + 20.0 - 20.0*std::exp(x: p1)-std::exp(x: p2);
67 }
68 return y;
69}
70
71Real sphere(const Array& x) {
72 //Minimum is found at 0
73 return DotProduct(v1: x, v2: x);
74}
75
76Array sphereValues(const Array& x) {
77 Array y(x.size());
78 for (Size i = 0; i < x.size(); i++) {
79 y[i] = x[i]*x[i];
80 }
81 return y;
82}
83
84Real rosenbrock(const Array& x) {
85 //Minimum is found at f(1, 1, ...)
86 QL_REQUIRE(x.size() > 1, "Input size needs to be higher than 1");
87 Real result = 0.0;
88 for (Size i = 0; i < x.size() - 1; i++) {
89 Real temp = (x[i + 1] - x[i] * x[i]);
90 result += (x[i] - 1.0)*(x[i] - 1.0) + 100.0*temp*temp;
91 }
92 return result;
93}
94
95Real easom(const Array& x) {
96 //Minimum is found at f(\pi, \pi, ...)
97 Real p1 = 1.0, p2 = 0.0;
98 for (Real i : x) {
99 p1 *= std::cos(x: i);
100 p2 += (i - M_PI) * (i - M_PI);
101 }
102 return -p1*std::exp(x: -p2);
103}
104
105Array easomValues(const Array& x) {
106 Array y(x.size());
107 for (Size i = 0; i < x.size(); i++) {
108 Real p1 = std::cos(x: x[i]);
109 Real p2 = (x[i] - M_PI)*(x[i] - M_PI);
110 y[i] = -p1*std::exp(x: -p2);
111 }
112 return y;
113}
114
115Real eggholder(const Array& x) {
116 //Minimum is found at f(512, 404.2319)
117 QL_REQUIRE(x.size() == 2, "Input size needs to be equal to 2");
118 Real p = (x[1] + 47.0);
119 return -p*std::sin(x: std::sqrt(x: std::abs(x: 0.5*x[0] + p))) -
120 x[0] * std::sin(x: std::sqrt(x: std::abs(x: x[0] - p)));
121}
122
123Real printFunction(Problem& p, const Array& x) {
124 std::cout << " f(" << x[0];
125 for (Size i = 1; i < x.size(); i++) {
126 std::cout << ", " << x[i];
127 }
128 Real val = p.value(x);
129 std::cout << ") = " << val << std::endl;
130 return val;
131}
132
133class TestFunction : public CostFunction {
134public:
135 typedef ext::function<Real(const Array&)> RealFunc;
136 typedef ext::function<Array(const Array&)> ArrayFunc;
137 explicit TestFunction(RealFunc f, ArrayFunc fs = ArrayFunc())
138 : f_(std::move(f)), fs_(std::move(fs)) {}
139 explicit TestFunction(Real (*f)(const Array&), Array (*fs)(const Array&) = nullptr)
140 : f_(f), fs_(fs) {}
141 ~TestFunction() override = default;
142 Real value(const Array& x) const override { return f_(x); }
143 Array values(const Array& x) const override {
144 if(!fs_)
145 throw std::runtime_error("Invalid function");
146 return fs_(x);
147 }
148
149private:
150 RealFunc f_;
151 ArrayFunc fs_;
152};
153
154int test(OptimizationMethod& method, CostFunction& f, const EndCriteria& endCriteria,
155 const Array& start, const Constraint& constraint = Constraint(),
156 const Array& optimum = Array()) {
157 QL_REQUIRE(!start.empty(), "Input size needs to be at least 1");
158 std::cout << "Starting point: ";
159 Constraint c;
160 if (!constraint.empty())
161 c = constraint;
162 Problem p(f, c, start);
163 printFunction(p, x: start);
164 method.minimize(P&: p, endCriteria);
165 std::cout << "End point: ";
166 Real val = printFunction(p, x: p.currentValue());
167 if(!optimum.empty())
168 {
169 std::cout << "Global optimum: ";
170 Real optimVal = printFunction(p, x: optimum);
171 if(std::abs(x: optimVal) < 1e-13)
172 return static_cast<int>(std::abs(x: val - optimVal) < 1e-6);
173 else
174 return static_cast<int>(std::abs(x: (val - optimVal) / optimVal) < 1e-6);
175 }
176 return 1;
177}
178
179void testFirefly() {
180 /*
181 The Eggholder function is only in 2 dimensions, it has a multitude
182 * of local minima, and they are not symmetric necessarily
183 */
184 Size n = 2;
185 NonhomogeneousBoundaryConstraint constraint(Array(n, -512.0), Array(n, 512.0));
186 Array x(n, 0.0);
187 Array optimum(n);
188 optimum[0] = 512.0;
189 optimum[1] = 404.2319;
190 Size agents = 150;
191 Real vola = 1.5;
192 Real intense = 1.0;
193 auto intensity = ext::make_shared<ExponentialIntensity>(args: 10.0, args: 1e-8, args&: intense);
194 auto randomWalk = ext::make_shared<LevyFlightWalk>(args&: vola, args: 0.5, args: 1.0, args&: seed);
195 std::cout << "Function eggholder, Agents: " << agents
196 << ", Vola: " << vola << ", Intensity: " << intense << std::endl;
197 TestFunction f(eggholder);
198 FireflyAlgorithm fa(agents, intensity, randomWalk, 40);
199 EndCriteria ec(5000, 1000, 1.0e-8, 1.0e-8, 1.0e-8);
200 test(method&: fa, f, endCriteria: ec, start: x, constraint, optimum);
201 std::cout << "================================================================" << std::endl;
202}
203
204void testSimulatedAnnealing(Size dimension, Size maxSteps, Size staticSteps){
205
206 /*The ackley function has a large amount of local minima, but the
207 structure is symmetric, so if one could simply just ignore the
208 walls separating the local minima, it would look like almost
209 like a parabola
210
211 Andres Hernandez: I could not find a configuration that was able
212 to fix the problem
213 */
214
215 //global minimum is at 0.0
216 TestFunction f(ackley, ackleyValues);
217
218 //Starting point
219 Array x(dimension, 1.5);
220 Array optimum(dimension, 0.0);
221
222 //Constraint for local optimizer
223 Array lower(dimension, -5.0);
224 Array upper(dimension, 5.0);
225 NonhomogeneousBoundaryConstraint constraint(lower, upper);
226
227 Real lambda = 0.1;
228 Real temperature = 350;
229 Real epsilon = 0.99;
230 Size ms = 1000;
231 std::cout << "Function ackley, Lambda: " << lambda
232 << ", Temperature: " << temperature
233 << ", Epsilon: " << epsilon
234 << ", Iterations: " << ms
235 << std::endl;
236
237 MersenneTwisterUniformRng rng(seed);
238 SimulatedAnnealing<MersenneTwisterUniformRng> sa(lambda, temperature, epsilon, ms, rng);
239 EndCriteria ec(maxSteps, staticSteps, 1.0e-8, 1.0e-8, 1.0e-8);
240 test(method&: sa, f, endCriteria: ec, start: x, constraint, optimum);
241 std::cout << "================================================================" << std::endl;
242}
243
244void testGaussianSA(Size dimension,
245 Size maxSteps,
246 Size staticSteps,
247 Real initialTemp,
248 Real finalTemp,
249 GaussianSimulatedAnnealing::ResetScheme resetScheme =
250 GaussianSimulatedAnnealing::ResetToBestPoint,
251 Size resetSteps = 150,
252 GaussianSimulatedAnnealing::LocalOptimizeScheme optimizeScheme =
253 GaussianSimulatedAnnealing::EveryBestPoint,
254 const ext::shared_ptr<OptimizationMethod>& localOptimizer =
255 ext::make_shared<LevenbergMarquardt>()) {
256
257 /*The ackley function has a large amount of local minima, but the
258 * structure is symmetric, so if one could simply just ignore the
259 * walls separating the local minima, it would look like almost like
260 * a parabola*/
261
262 //global minimum is at 0.0
263 TestFunction f(ackley, ackleyValues);
264
265 std::cout << "Function: ackley, Dimensions: " << dimension
266 << ", Initial temp:" << initialTemp
267 << ", Final temp:" << finalTemp
268 << ", Reset scheme:" << resetScheme
269 << ", Reset steps:" << resetSteps
270 << std::endl;
271 //Starting point
272 Array x(dimension, 1.5);
273 Array optimum(dimension, 0.0);
274
275 //Constraint for local optimizer
276 Array lower(dimension, -5.0);
277 Array upper(dimension, 5.0);
278 NonhomogeneousBoundaryConstraint constraint(lower, upper);
279
280 //Simulated annealing setup
281 SamplerGaussian sampler(seed);
282 ProbabilityBoltzmannDownhill probability(seed);
283 TemperatureExponential temperature(initialTemp, dimension);
284 GaussianSimulatedAnnealing sa(sampler, probability, temperature, ReannealingTrivial(),
285 initialTemp, finalTemp, 50, resetScheme,
286 resetSteps, localOptimizer,
287 optimizeScheme);
288
289 EndCriteria ec(maxSteps, staticSteps, 1.0e-8, 1.0e-8, 1.0e-8);
290 test(method&: sa, f, endCriteria: ec, start: x, constraint, optimum);
291 std::cout << "================================================================" << std::endl;
292}
293
294void testPSO(Size n){
295 /*The Rosenbrock function has a global minima at (1.0, ...) and a local minima at (-1.0, 1.0, ...)
296 The difficulty lies in the weird shape of the function*/
297 NonhomogeneousBoundaryConstraint constraint(Array(n, -1.0), Array(n, 4.0));
298 Array x(n, 0.0);
299 Array optimum(n, 1.0);
300 Size agents = 100;
301 Size kneighbor = 25;
302 Size threshold = 500;
303 std::cout << "Function: rosenbrock, Dimensions: " << n
304 << ", Agents: " << agents << ", K-neighbors: " << kneighbor
305 << ", Threshold: " << threshold << std::endl;
306 auto topology = ext::make_shared<KNeighbors>(args&: kneighbor);
307 auto inertia = ext::make_shared<LevyFlightInertia>(args: 1.5, args&: threshold, args&: seed);
308 TestFunction f(rosenbrock);
309 ParticleSwarmOptimization pso(agents, topology, inertia, 2.05, 2.05, seed);
310 EndCriteria ec(10000, 1000, 1.0e-8, 1.0e-8, 1.0e-8);
311 test(method&: pso, f, endCriteria: ec, start: x, constraint, optimum);
312 std::cout << "================================================================" << std::endl;
313}
314
315void testDifferentialEvolution(Size n, Size agents){
316 /*The Rosenbrock function has a global minima at (1.0, ...) and a local minima at (-1.0, 1.0, ...)
317 The difficulty lies in the weird shape of the function*/
318 NonhomogeneousBoundaryConstraint constraint(Array(n, -4.0), Array(n, 4.0));
319 Array x(n, 0.0);
320 Array optimum(n, 1.0);
321
322 TestFunction f(rosenbrock);
323
324 Real probability = 0.3;
325 Real stepsizeWeight = 0.6;
326 DifferentialEvolution::Strategy strategy = DifferentialEvolution::BestMemberWithJitter;
327
328 std::cout << "Function: rosenbrock, Dimensions: " << n << ", Agents: " << agents
329 << ", Probability: " << probability
330 << ", StepsizeWeight: " << stepsizeWeight
331 << ", Strategy: BestMemberWithJitter" << std::endl;
332 DifferentialEvolution::Configuration config;
333 config.withBounds(b: true)
334 .withCrossoverProbability(p: probability)
335 .withPopulationMembers(n: agents)
336 .withStepsizeWeight(w: stepsizeWeight)
337 .withStrategy(s: strategy)
338 .withSeed(s: seed);
339
340 DifferentialEvolution de(config);
341 EndCriteria ec(5000, 1000, 1.0e-8, 1.0e-8, 1.0e-8);
342 test(method&: de, f, endCriteria: ec, start: x, constraint, optimum);
343 std::cout << "================================================================" << std::endl;
344}
345
346
347int main(int, char* []) {
348
349 try {
350 std::cout << std::endl;
351
352 std::cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
353 std::cout << "Firefly Algorithm Test" << std::endl;
354 std::cout << "----------------------------------------------------------------" << std::endl;
355 testFirefly();
356
357
358 std::cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
359 std::cout << "Hybrid Simulated Annealing Test" << std::endl;
360 std::cout << "----------------------------------------------------------------" << std::endl;
361 testGaussianSA(dimension: 3, maxSteps: 500, staticSteps: 200, initialTemp: 100.0, finalTemp: 0.1, resetScheme: GaussianSimulatedAnnealing::ResetToBestPoint, resetSteps: 150, optimizeScheme: GaussianSimulatedAnnealing::EveryNewPoint);
362 testGaussianSA(dimension: 10, maxSteps: 500, staticSteps: 200, initialTemp: 100.0, finalTemp: 0.1, resetScheme: GaussianSimulatedAnnealing::ResetToBestPoint, resetSteps: 150, optimizeScheme: GaussianSimulatedAnnealing::EveryNewPoint);
363 testGaussianSA(dimension: 30, maxSteps: 500, staticSteps: 200, initialTemp: 100.0, finalTemp: 0.1, resetScheme: GaussianSimulatedAnnealing::ResetToBestPoint, resetSteps: 150, optimizeScheme: GaussianSimulatedAnnealing::EveryNewPoint);
364
365
366 std::cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
367 std::cout << "Particle Swarm Optimization Test" << std::endl;
368 std::cout << "----------------------------------------------------------------" << std::endl;
369 testPSO(n: 3);
370 testPSO(n: 10);
371 testPSO(n: 30);
372
373
374 std::cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
375 std::cout << "Simulated Annealing Test" << std::endl;
376 std::cout << "----------------------------------------------------------------" << std::endl;
377 testSimulatedAnnealing(dimension: 3, maxSteps: 10000, staticSteps: 4000);
378 testSimulatedAnnealing(dimension: 10, maxSteps: 10000, staticSteps: 4000);
379 testSimulatedAnnealing(dimension: 30, maxSteps: 10000, staticSteps: 4000);
380
381
382 std::cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
383 std::cout << "Differential Evolution Test" << std::endl;
384 std::cout << "----------------------------------------------------------------" << std::endl;
385 testDifferentialEvolution(n: 3, agents: 50);
386 testDifferentialEvolution(n: 10, agents: 150);
387 testDifferentialEvolution(n: 30, agents: 450);
388
389 return 0;
390 } catch (std::exception& e) {
391 std::cerr << e.what() << std::endl;
392 return 1;
393 } catch (...) {
394 std::cerr << "unknown error" << std::endl;
395 return 1;
396 }
397}
398

source code of quantlib/Examples/GlobalOptimizer/GlobalOptimizer.cpp

Morty Proxy This is a proxified and sanitized view of the page, visit original site.