1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2001, 2002, 2003 Sadruddin Rejeb
5 Copyright (C) 2006 Ferdinando Ametrano
6 Copyright (C) 2007 Marco Bianchetti
7 Copyright (C) 2007 François du Vignaud
8
9 This file is part of QuantLib, a free-software/open-source library
10 for financial quantitative analysts and developers - http://quantlib.org/
11
12 QuantLib is free software: you can redistribute it and/or modify it
13 under the terms of the QuantLib license. You should have received a
14 copy of the license along with this program; if not, please email
15 <quantlib-dev@lists.sf.net>. The license is also available online at
16 <http://quantlib.org/license.shtml>.
17
18 This program is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20 FOR A PARTICULAR PURPOSE. See the license for more details.
21*/
22
23/* The implementation of the algorithm was highly inspired by
24 * "Numerical Recipes in C", 2nd edition, Press, Teukolsky, Vetterling,
25 * Flannery, chapter 10.
26 * Modified may 2007: end criteria set on x instead on fx,
27 * inspired by bad behaviour found with test function fx=x*x+x+1,
28 * xStart = -100, lambda = 1.0, ftol = 1.e-16
29 * (it reports x=0 as the minimum!)
30 * and by GSL implementation, v. 1.9 (http://www.gnu.org/software/gsl/)
31 */
32
33#include <ql/math/optimization/simplex.hpp>
34#include <ql/math/optimization/constraint.hpp>
35
36#if !defined(__GNUC__) || __GNUC__ > 3 || __GNUC_MINOR__ > 4
37#define QL_ARRAY_EXPRESSIONS
38#endif
39
40namespace QuantLib {
41
42 namespace {
43 // Computes the size of the simplex
44 Real computeSimplexSize (const std::vector<Array>& vertices) {
45 Array center(vertices.front().size(),0);
46 for (const auto& vertice : vertices)
47 center += vertice;
48 center *=1/Real(vertices.size());
49 Real result = 0;
50 for (const auto& vertice : vertices) {
51 Array temp = vertice - center;
52 result += Norm2(v: temp);
53 }
54 return result/Real(vertices.size());
55 }
56 }
57
58 Real Simplex::extrapolate(Problem& P,
59 Size iHighest,
60 Real &factor) const {
61
62 Array pTry;
63 do {
64 Size dimensions = values_.size() - 1;
65 Real factor1 = (1.0 - factor)/dimensions;
66 Real factor2 = factor1 - factor;
67 #if defined(QL_ARRAY_EXPRESSIONS)
68 pTry = sum_*factor1 - vertices_[iHighest]*factor2;
69 #else
70 // composite expressions fail to compile with gcc 3.4 on windows
71 pTry = sum_*factor1;
72 pTry -= vertices_[iHighest]*factor2;
73 #endif
74 factor *= 0.5;
75 } while (!P.constraint().test(p: pTry) && std::fabs(x: factor) > QL_EPSILON);
76 if (std::fabs(x: factor) <= QL_EPSILON) {
77 return values_[iHighest];
78 }
79 factor *= 2.0;
80 Real vTry = P.value(x: pTry);
81 if (vTry < values_[iHighest]) {
82 values_[iHighest] = vTry;
83 #if defined(QL_ARRAY_EXPRESSIONS)
84 sum_ += pTry - vertices_[iHighest];
85 #else
86 sum_ += pTry;
87 sum_ -= vertices_[iHighest];
88 #endif
89 vertices_[iHighest] = pTry;
90 }
91 return vTry;
92
93 }
94
95
96 EndCriteria::Type Simplex::minimize(Problem& P,
97 const EndCriteria& endCriteria) {
98 // set up of the problem
99 //Real ftol = endCriteria.functionEpsilon(); // end criteria on f(x) (see Numerical Recipes in C++, p.410)
100 Real xtol = endCriteria.rootEpsilon(); // end criteria on x (see GSL v. 1.9, http://www.gnu.org/software/gsl/)
101 Size maxStationaryStateIterations_
102 = endCriteria.maxStationaryStateIterations();
103 EndCriteria::Type ecType = EndCriteria::None;
104 P.reset();
105
106 Array x_ = P.currentValue();
107 if (!P.constraint().test(p: x_))
108 QL_FAIL("Initial guess " << x_ << " is not in the feasible region.");
109
110 Integer iterationNumber_=0;
111
112 // Initialize vertices of the simplex
113 Size n = x_.size();
114 vertices_ = std::vector<Array>(n+1, x_);
115 for (Size i=0; i<n; ++i) {
116 Array direction(n, 0.0);
117 direction[i] = 1.0;
118 P.constraint().update(p&: vertices_[i+1], direction, beta: lambda_);
119 }
120 // Initialize function values at the vertices of the simplex
121 values_ = Array(n+1, 0.0);
122 for (Size i=0; i<=n; ++i)
123 values_[i] = P.value(x: vertices_[i]);
124 // Loop looking for minimum
125 do {
126 sum_ = Array(n, 0.0);
127 Size i;
128 for (i=0; i<=n; i++)
129 sum_ += vertices_[i];
130 // Determine the best (iLowest), worst (iHighest)
131 // and 2nd worst (iNextHighest) vertices
132 Size iLowest = 0;
133 Size iHighest, iNextHighest;
134 if (values_[0]<values_[1]) {
135 iHighest = 1;
136 iNextHighest = 0;
137 } else {
138 iHighest = 0;
139 iNextHighest = 1;
140 }
141 for (i=1;i<=n; i++) {
142 if (values_[i]>values_[iHighest]) {
143 iNextHighest = iHighest;
144 iHighest = i;
145 } else {
146 if ((values_[i]>values_[iNextHighest]) && i!=iHighest)
147 iNextHighest = i;
148 }
149 if (values_[i]<values_[iLowest])
150 iLowest = i;
151 }
152 // Now compute accuracy, update iteration number and check end criteria
153 //// Numerical Recipes exit strategy on fx (see NR in C++, p.410)
154 //Real low = values_[iLowest];
155 //Real high = values_[iHighest];
156 //Real rtol = 2.0*std::fabs(high - low)/
157 // (std::fabs(high) + std::fabs(low) + QL_EPSILON);
158 //++iterationNumber_;
159 //if (rtol < ftol ||
160 // endCriteria.checkMaxIterations(iterationNumber_, ecType)) {
161 // GSL exit strategy on x (see GSL v. 1.9, http://www.gnu.org/software/gsl
162 Real simplexSize = computeSimplexSize(vertices: vertices_);
163 ++iterationNumber_;
164 if (simplexSize < xtol ||
165 endCriteria.checkMaxIterations(iteration: iterationNumber_, ecType)) {
166 endCriteria.checkStationaryPoint(xOld: 0.0, xNew: 0.0,
167 statStateIterations&: maxStationaryStateIterations_, ecType);
168 endCriteria.checkMaxIterations(iteration: iterationNumber_, ecType);
169 x_ = vertices_[iLowest];
170 Real low = values_[iLowest];
171 P.setFunctionValue(low);
172 P.setCurrentValue(x_);
173 return ecType;
174 }
175 // If end criteria is not met, continue
176 Real factor = -1.0;
177 Real vTry = extrapolate(P, iHighest, factor);
178 if ((vTry <= values_[iLowest]) && (factor == -1.0)) {
179 factor = 2.0;
180 extrapolate(P, iHighest, factor);
181 } else if (std::fabs(x: factor) > QL_EPSILON) {
182 if (vTry >= values_[iNextHighest]) {
183 Real vSave = values_[iHighest];
184 factor = 0.5;
185 vTry = extrapolate(P, iHighest, factor);
186 if (vTry >= vSave && std::fabs(x: factor) > QL_EPSILON) {
187 for (Size i=0; i<=n; i++) {
188 if (i!=iLowest) {
189 #if defined(QL_ARRAY_EXPRESSIONS)
190 vertices_[i] =
191 0.5*(vertices_[i] + vertices_[iLowest]);
192 #else
193 vertices_[i] += vertices_[iLowest];
194 vertices_[i] *= 0.5;
195 #endif
196 values_[i] = P.value(x: vertices_[i]);
197 }
198 }
199 }
200 }
201 }
202 // If can't extrapolate given the constraints, exit
203 if (std::fabs(x: factor) <= QL_EPSILON) {
204 x_ = vertices_[iLowest];
205 Real low = values_[iLowest];
206 P.setFunctionValue(low);
207 P.setCurrentValue(x_);
208 return EndCriteria::StationaryFunctionValue;
209 }
210 } while (true);
211 QL_FAIL("optimization failed: unexpected behaviour");
212 }
213}
214

source code of quantlib/ql/math/optimization/simplex.cpp