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