1 |
tim |
1002 |
#include "Minimizer1D.hpp"
|
2 |
tim |
1005 |
#include "math.h"
|
3 |
|
|
|
4 |
|
|
GoldenSectionMinimizer::GoldenSectionMinimizer(NLModel* nlp)
|
5 |
|
|
:Minimizer1D(nlp){
|
6 |
|
|
setName("GoldenSection");
|
7 |
|
|
}
|
8 |
tim |
1002 |
|
9 |
tim |
1005 |
int GoldenSectionMinimizer::checkConvergence(){
|
10 |
|
|
|
11 |
tim |
1002 |
if ((rightVar - leftVar) < stepTol)
|
12 |
tim |
1010 |
return 1;
|
13 |
tim |
1002 |
else
|
14 |
|
|
return -1;
|
15 |
|
|
}
|
16 |
|
|
void GoldenSectionMinimizer::minimize(){
|
17 |
|
|
vector<double> tempX;
|
18 |
|
|
vector <double> currentX;
|
19 |
|
|
|
20 |
|
|
const double goldenRatio = 0.618034;
|
21 |
|
|
|
22 |
tim |
1015 |
tempX = currentX = model->getX();
|
23 |
tim |
1002 |
|
24 |
|
|
alpha = leftVar + (1 - goldenRatio) * (rightVar - leftVar);
|
25 |
|
|
beta = leftVar + goldenRatio * (rightVar - leftVar);
|
26 |
|
|
|
27 |
tim |
1015 |
for (int i = 0; i < tempX.size(); i ++)
|
28 |
|
|
tempX[i] = currentX[i] + direction[i] * alpha;
|
29 |
|
|
|
30 |
tim |
1002 |
fAlpha = model->calcF(tempX);
|
31 |
|
|
|
32 |
tim |
1015 |
for (int i = 0; i < tempX.size(); i ++)
|
33 |
|
|
tempX[i] = currentX[i] + direction[i] * beta;
|
34 |
|
|
|
35 |
tim |
1002 |
fBeta = model->calcF(tempX);
|
36 |
|
|
|
37 |
|
|
for(currentIter = 0; currentIter < maxIteration; currentIter++){
|
38 |
|
|
|
39 |
|
|
if (checkConvergence() > 0){
|
40 |
|
|
minStatus = MINSTATUS_CONVERGE;
|
41 |
|
|
return;
|
42 |
|
|
}
|
43 |
|
|
|
44 |
|
|
if (fAlpha > fBeta){
|
45 |
|
|
leftVar = alpha;
|
46 |
|
|
alpha = beta;
|
47 |
|
|
beta = leftVar + goldenRatio * (rightVar - leftVar);
|
48 |
|
|
|
49 |
tim |
1015 |
for (int i = 0; i < tempX.size(); i ++)
|
50 |
|
|
tempX[i] = currentX[i] + direction[i] * beta;
|
51 |
tim |
1002 |
|
52 |
|
|
prevMinVar = minVar;
|
53 |
|
|
fPrevMinVar = fMinVar;
|
54 |
|
|
|
55 |
|
|
minVar = beta;
|
56 |
|
|
fMinVar = model->calcF(tempX);
|
57 |
|
|
|
58 |
|
|
}
|
59 |
|
|
else{
|
60 |
|
|
rightVar = beta;
|
61 |
|
|
beta = alpha;
|
62 |
|
|
alpha = leftVar + (1 - goldenRatio) * (rightVar - leftVar);
|
63 |
|
|
|
64 |
tim |
1015 |
for (int i = 0; i < tempX.size(); i ++)
|
65 |
|
|
tempX[i] = currentX[i] + direction[i] * alpha;
|
66 |
tim |
1002 |
|
67 |
|
|
prevMinVar = minVar;
|
68 |
|
|
fPrevMinVar = fMinVar;
|
69 |
|
|
|
70 |
|
|
minVar = alpha;
|
71 |
|
|
fMinVar = model->calcF(tempX);
|
72 |
|
|
}
|
73 |
|
|
|
74 |
|
|
}
|
75 |
|
|
|
76 |
|
|
minStatus = MINSTATUS_MAXITER;
|
77 |
|
|
|
78 |
|
|
}
|
79 |
|
|
|
80 |
tim |
1005 |
/**
|
81 |
|
|
* Brent's method is a root-finding algorithm which combines root bracketing, interval bisection,
|
82 |
|
|
* and inverse quadratic interpolation.
|
83 |
tim |
1002 |
*/
|
84 |
tim |
1005 |
BrentMinimizer::BrentMinimizer(NLModel* nlp)
|
85 |
|
|
:Minimizer1D(nlp){
|
86 |
|
|
setName("Brent");
|
87 |
|
|
}
|
88 |
tim |
1002 |
|
89 |
|
|
void BrentMinimizer::minimize(){
|
90 |
|
|
|
91 |
tim |
1005 |
double fu, fv, fw;
|
92 |
|
|
double p, q, r;
|
93 |
|
|
double u, v, w;
|
94 |
|
|
double d;
|
95 |
|
|
double e;
|
96 |
|
|
double etemp;
|
97 |
|
|
double stepTol2;
|
98 |
tim |
1010 |
double fLeftVar, fRightVar;
|
99 |
tim |
1005 |
const double goldenRatio = 0.3819660;
|
100 |
|
|
vector<double> tempX, currentX;
|
101 |
|
|
|
102 |
|
|
stepTol2 = 2 * stepTol;
|
103 |
|
|
e = 0;
|
104 |
|
|
d = 0;
|
105 |
|
|
|
106 |
|
|
currentX = tempX = model->getX();
|
107 |
|
|
|
108 |
tim |
1015 |
for (int i = 0; i < tempX.size(); i ++)
|
109 |
|
|
tempX[i] = currentX[i] + direction[i] * leftVar;
|
110 |
|
|
|
111 |
tim |
1010 |
fLeftVar = model->calcF(tempX);
|
112 |
tim |
1015 |
|
113 |
|
|
for (int i = 0; i < tempX.size(); i ++)
|
114 |
|
|
tempX[i] = currentX[i] + direction[i] * rightVar;
|
115 |
tim |
1005 |
|
116 |
tim |
1010 |
fRightVar = model->calcF(tempX);
|
117 |
tim |
1005 |
|
118 |
tim |
1010 |
if(fRightVar < fLeftVar) {
|
119 |
|
|
prevMinVar = rightVar;
|
120 |
|
|
fPrevMinVar = fRightVar;
|
121 |
tim |
1005 |
v = leftVar;
|
122 |
|
|
fv = fLeftVar;
|
123 |
|
|
}
|
124 |
|
|
else {
|
125 |
|
|
prevMinVar = leftVar;
|
126 |
tim |
1010 |
fPrevMinVar = fLeftVar;
|
127 |
tim |
1005 |
v = rightVar;
|
128 |
tim |
1010 |
fv = fRightVar;
|
129 |
tim |
1005 |
}
|
130 |
tim |
1010 |
|
131 |
|
|
midVar = leftVar + rightVar;
|
132 |
tim |
1005 |
|
133 |
tim |
1002 |
for(currentIter = 0; currentIter < maxIteration; currentIter){
|
134 |
|
|
|
135 |
tim |
1005 |
// a trial parabolic fit
|
136 |
|
|
if (fabs(e) > stepTol){
|
137 |
tim |
1002 |
|
138 |
tim |
1005 |
r = (minVar - prevMinVar) * (fMinVar - fv);
|
139 |
|
|
q = (minVar - v) * (fMinVar - fPrevMinVar);
|
140 |
|
|
p = (minVar - v) *q -(minVar - prevMinVar)*r;
|
141 |
|
|
q = 2.0 *(q-r);
|
142 |
tim |
1002 |
|
143 |
tim |
1005 |
if (q > 0.0)
|
144 |
|
|
p = -p;
|
145 |
tim |
1002 |
|
146 |
tim |
1005 |
q = fabs(q);
|
147 |
|
|
|
148 |
|
|
etemp = e;
|
149 |
|
|
e = d;
|
150 |
|
|
|
151 |
tim |
1010 |
if(fabs(p) >= fabs(0.5*q*etemp) || p <= q*(leftVar - minVar) || p >= q*(rightVar - minVar)){
|
152 |
tim |
1005 |
e = minVar >= midVar ? leftVar - minVar : rightVar - minVar;
|
153 |
|
|
d = goldenRatio * e;
|
154 |
|
|
}
|
155 |
|
|
else{
|
156 |
|
|
d = p/q;
|
157 |
|
|
u = minVar + d;
|
158 |
|
|
if ( u - leftVar < stepTol2 || rightVar - u < stepTol2)
|
159 |
|
|
d = midVar > minVar ? stepTol : - stepTol;
|
160 |
|
|
}
|
161 |
|
|
}
|
162 |
|
|
//golden section
|
163 |
|
|
else{
|
164 |
|
|
e = minVar >=midVar? leftVar - minVar : rightVar - minVar;
|
165 |
|
|
d =goldenRatio * e;
|
166 |
|
|
}
|
167 |
|
|
|
168 |
|
|
u = fabs(d) >= stepTol ? minVar + d : minVar + copysign(d, stepTol);
|
169 |
|
|
|
170 |
tim |
1015 |
for (int i = 0; i < tempX.size(); i ++)
|
171 |
|
|
tempX[i] = currentX[i] + direction[i] * u;
|
172 |
|
|
|
173 |
tim |
1005 |
fu = model->calcF();
|
174 |
|
|
|
175 |
|
|
if(fu <= fMinVar){
|
176 |
|
|
|
177 |
|
|
if(u >= minVar)
|
178 |
|
|
leftVar = minVar;
|
179 |
|
|
else
|
180 |
|
|
rightVar = minVar;
|
181 |
|
|
|
182 |
|
|
v = prevMinVar;
|
183 |
|
|
fv = fPrevMinVar;
|
184 |
|
|
prevMinVar = minVar;
|
185 |
|
|
fPrevMinVar = fMinVar;
|
186 |
|
|
minVar = u;
|
187 |
|
|
fMinVar = fu;
|
188 |
|
|
|
189 |
|
|
}
|
190 |
|
|
else{
|
191 |
|
|
if (u < minVar) leftVar = u;
|
192 |
|
|
else rightVar= u;
|
193 |
|
|
if(fu <= fPrevMinVar || prevMinVar == minVar) {
|
194 |
|
|
v = prevMinVar;
|
195 |
|
|
fv = fPrevMinVar;
|
196 |
|
|
prevMinVar = u;
|
197 |
|
|
fPrevMinVar = fu;
|
198 |
|
|
}
|
199 |
|
|
else if ( fu <= fv || v == minVar || v == prevMinVar ) {
|
200 |
|
|
v = u;
|
201 |
|
|
fv = fu;
|
202 |
|
|
}
|
203 |
|
|
}
|
204 |
|
|
|
205 |
|
|
midVar = leftVar + rightVar;
|
206 |
|
|
|
207 |
|
|
if (checkConvergence() > 0){
|
208 |
|
|
minStatus = MINSTATUS_CONVERGE;
|
209 |
|
|
return;
|
210 |
|
|
}
|
211 |
|
|
|
212 |
tim |
1002 |
}
|
213 |
|
|
|
214 |
|
|
|
215 |
|
|
minStatus = MINSTATUS_MAXITER;
|
216 |
tim |
1010 |
return;
|
217 |
tim |
1002 |
}
|
218 |
tim |
1005 |
|
219 |
tim |
1010 |
int BrentMinimizer::checkConvergence(){
|
220 |
tim |
1005 |
|
221 |
|
|
if (fabs(minVar - midVar) < stepTol)
|
222 |
|
|
return 1;
|
223 |
|
|
else
|
224 |
|
|
return -1;
|
225 |
|
|
}
|