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