| 1 |
|
#include "Minimizer1D.hpp" |
| 2 |
+ |
#include "math.h" |
| 3 |
+ |
|
| 4 |
+ |
//----------------------------------------------------------------------------// |
| 5 |
|
void Minimizer1D::Minimize(vector<double>& direction), double left, double right); { |
| 6 |
|
setDirection(direction); |
| 7 |
|
setRange(left,right); |
| 8 |
|
minimize(); |
| 9 |
|
} |
| 10 |
|
|
| 11 |
< |
int Minimizer1D::checkConvergence(){ |
| 11 |
> |
//----------------------------------------------------------------------------// |
| 12 |
> |
GoldenSectionMinimizer::GoldenSectionMinimizer(NLModel* nlp) |
| 13 |
> |
:Minimizer1D(nlp){ |
| 14 |
> |
setName("GoldenSection"); |
| 15 |
> |
} |
| 16 |
|
|
| 17 |
+ |
int GoldenSectionMinimizer::checkConvergence(){ |
| 18 |
+ |
|
| 19 |
|
if ((rightVar - leftVar) < stepTol) |
| 20 |
< |
return |
| 20 |
> |
return 1 |
| 21 |
|
else |
| 22 |
|
return -1; |
| 23 |
|
} |
| 15 |
– |
|
| 24 |
|
void GoldenSectionMinimizer::minimize(){ |
| 25 |
|
vector<double> tempX; |
| 26 |
|
vector <double> currentX; |
| 79 |
|
|
| 80 |
|
} |
| 81 |
|
|
| 82 |
< |
/* |
| 83 |
< |
* |
| 82 |
> |
/** |
| 83 |
> |
* Brent's method is a root-finding algorithm which combines root bracketing, interval bisection, |
| 84 |
> |
* and inverse quadratic interpolation. |
| 85 |
|
*/ |
| 86 |
+ |
BrentMinimizer::BrentMinimizer(NLModel* nlp) |
| 87 |
+ |
:Minimizer1D(nlp){ |
| 88 |
+ |
setName("Brent"); |
| 89 |
+ |
} |
| 90 |
|
|
| 91 |
|
void BrentMinimizer::minimize(){ |
| 92 |
|
|
| 93 |
+ |
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 |
|
for(currentIter = 0; currentIter < maxIteration; currentIter){ |
| 130 |
|
|
| 131 |
+ |
// a trial parabolic fit |
| 132 |
+ |
if (fabs(e) > stepTol){ |
| 133 |
|
|
| 134 |
+ |
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 |
|
|
| 139 |
+ |
if (q > 0.0) |
| 140 |
+ |
p = -p; |
| 141 |
|
|
| 142 |
+ |
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 |
|
} |
| 207 |
|
|
| 208 |
|
|
| 209 |
|
minStatus = MINSTATUS_MAXITER; |
| 210 |
|
return; |
| 211 |
+ |
|
| 212 |
+ |
//----------------------------------------------------------------------------// |
| 213 |
+ |
|
| 214 |
|
} |
| 215 |
+ |
|
| 216 |
+ |
BrentMinimizer::checkConvergence(){ |
| 217 |
+ |
|
| 218 |
+ |
if (fabs(minVar - midVar) < stepTol) |
| 219 |
+ |
return 1; |
| 220 |
+ |
else |
| 221 |
+ |
return -1; |
| 222 |
+ |
} |