ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/Minimizer1D.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/Minimizer1D.cpp (file contents):
Revision 1002 by tim, Mon Feb 2 20:29:41 2004 UTC vs.
Revision 1031 by tim, Fri Feb 6 18:58:06 2004 UTC

# Line 1 | Line 1
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   }
# Line 19 | Line 20 | void GoldenSectionMinimizer::minimize(){
20  
21    const double goldenRatio = 0.618034;
22    
23 <  currentX =  model->getX();
23 >  tempX = currentX =  model->getX();
24  
25    alpha = leftVar + (1 - goldenRatio) * (rightVar  - leftVar);
26    beta = leftVar + goldenRatio * (rightVar - leftVar);
27  
28 <  tempX = currentX + direction * alpha;
28 >  for (int i = 0; i < tempX.size(); i ++)
29 >    tempX[i] = currentX[i] + direction[i] * alpha;
30 >
31    fAlpha = model->calcF(tempX);
32  
33 <  tempX = currentX + direction * beta;
33 >  for (int i = 0; i < tempX.size(); i ++)
34 >    tempX[i] = currentX[i] + direction[i] * beta;
35 >
36    fBeta = model->calcF(tempX);
37  
38    for(currentIter = 0; currentIter < maxIteration; currentIter++){
# Line 40 | Line 45 | void GoldenSectionMinimizer::minimize(){
45      if (fAlpha > fBeta){
46        leftVar = alpha;
47        alpha = beta;
48 +      
49        beta =  leftVar + goldenRatio * (rightVar - leftVar);
50  
51 <      tempX = currentX + beta * direction;
52 <
53 <      prevMinVar = minVar;
54 <      fPrevMinVar = fMinVar;
51 >      for (int i = 0; i < tempX.size(); i ++)
52 >        tempX[i] = currentX[i] + direction[i] * beta;
53 >      fAlpha = fBeta;
54 >      fBeta  = model->calcF(tempX);
55        
56 +      prevMinVar = alpha;
57 +      fPrevMinVar = fAlpha;      
58        minVar = beta;
59 <      fMinVar = model->calcF(tempX);
52 <      
59 >      fMinVar = fBeta;      
60      }
61      else{
62        rightVar = beta;
63        beta = alpha;
64 +
65        alpha = leftVar + (1 - goldenRatio) * (rightVar  - leftVar);
66  
67 <      tempX = currentX + alpha * direction;
67 >      for (int i = 0; i < tempX.size(); i ++)
68 >        tempX[i] = currentX[i] + direction[i] * alpha;
69  
70 <      prevMinVar = minVar;
71 <      fPrevMinVar = fMinVar;
72 <      
70 >      fBeta = fAlpha;
71 >      fAlpha = model->calcF(tempX);
72 >
73 >      prevMinVar = beta;
74 >      fPrevMinVar = fBeta;
75 >
76        minVar = alpha;
77 <      fMinVar = model->calcF(tempX);
77 >      fMinVar = fAlpha;      
78      }
79      
80    }
# Line 71 | Line 83 | void GoldenSectionMinimizer::minimize(){
83      
84   }
85  
86 < /*
87 < *
86 > /**
87 > * Brent's method is a root-finding algorithm which combines root bracketing, interval bisection,
88 > * and inverse quadratic interpolation.
89   */
90 + BrentMinimizer::BrentMinimizer(NLModel* nlp)
91 +                   :Minimizer1D(nlp){
92 +  setName("Brent");
93 + }
94  
95 + void BrentMinimizer::minimize(vector<double>& direction, double left, double right){
96 +
97 +  //brent algorithm ascending order
98 +
99 +  if (left > right)
100 +      setRange(right, left);
101 +  else    
102 +    setRange(left, right);
103 +  
104 +  setDirection(direction);
105 +  
106 +  minimize();
107 + }
108   void BrentMinimizer::minimize(){
109  
110 <  for(currentIter = 0; currentIter < maxIteration; currentIter){
110 >  double fu, fv, fw;
111 >  double p, q, r;
112 >  double u, v, w;
113 >  double d;
114 >  double e;
115 >  double etemp;
116 >  double stepTol2;
117 >  double fLeftVar, fRightVar;
118 >  const double goldenRatio = 0.3819660;
119 >  vector<double> tempX, currentX;
120 >  
121 >  stepTol2 = 2 * stepTol;
122  
123 +  e = 0;
124 +  d = 0;
125  
126 +  currentX = model->getX();
127 +  tempX.resize(currentX.size());
128  
129 +  
130  
131 +  
132 +  for (int i = 0; i < tempX.size(); i ++)
133 +    tempX[i] = currentX[i] + direction[i] * leftVar;
134 +  
135 +  fLeftVar = model->calcF(tempX);
136 +
137 +  for (int i = 0; i < tempX.size(); i ++)
138 +    tempX[i] = currentX[i] + direction[i] * rightVar;  
139 +  
140 +  fRightVar = model->calcF(tempX);
141 +
142 +  // find an interior point  left < interior < right which satisfy f(left) > f(interior) and f(right) > f(interior)
143 +  
144 +  bracket(minVar, fMinVar, leftVar, fLeftVar, rightVar, fRightVar);
145 +
146 +  if(fRightVar < fLeftVar) {
147 +    prevMinVar = rightVar;
148 +    fPrevMinVar = fRightVar;
149 +    v  = leftVar;
150 +    fv = fLeftVar;
151    }
152 +  else {
153 +    prevMinVar = leftVar;
154 +    fPrevMinVar = fLeftVar;
155 +    v  = rightVar;
156 +    fv = fRightVar;
157 +  }
158 +    
159 +  minVar =  rightVar+ goldenRatio * (rightVar - leftVar);
160  
161 +  for (int i = 0; i < tempX.size(); i ++)
162 +    tempX[i] = currentX[i] + direction[i] * minVar;
163  
164 +  fMinVar = model->calcF(tempX);
165 +  
166 +  prevMinVar = v = minVar;
167 +  fPrevMinVar = fv = fMinVar;
168 +  midVar = (leftVar + rightVar) / 2;
169 +  
170 +  for(currentIter = 0; currentIter < maxIteration; currentIter++){
171 +
172 +     //construct a trial parabolic fit
173 +     if (fabs(e) > stepTol){
174 +
175 +       r = (minVar - prevMinVar) * (fMinVar - fv);
176 +       q = (minVar - v) * (fMinVar - fPrevMinVar);
177 +       p = (minVar - v) *q -(minVar - prevMinVar)*r;
178 +       q = 2.0 *(q-r);
179 +
180 +       if (q > 0.0)
181 +         p = -p;
182 +
183 +       q = fabs(q);
184 +
185 +       etemp = e;
186 +       e  = d;
187 +
188 +       if(fabs(p) >= fabs(0.5*q*etemp) || p <= q*(leftVar - minVar) || p >= q*(rightVar - minVar)){
189 +         //reject parabolic fit and use golden section step instead
190 +         e =  minVar >= midVar ? leftVar - minVar : rightVar - minVar;
191 +         d = goldenRatio * e;
192 +       }
193 +       else{
194 +
195 +        //take the parabolic step
196 +         d = p/q;
197 +         u = minVar + d;
198 +         if ( u - leftVar < stepTol2 || rightVar - u  < stepTol2)
199 +           d = midVar > minVar ? stepTol : - stepTol;
200 +       }
201 +      
202 +     }
203 +     else{
204 +       e = minVar >= midVar ? leftVar -minVar : rightVar-minVar;
205 +       d = goldenRatio * e;        
206 +     }
207 +
208 +     u = fabs(d) >= stepTol ? minVar + d : minVar + copysign(stepTol, d);
209 +
210 +     for (int i = 0; i < tempX.size(); i ++)
211 +       tempX[i] = currentX[i] + direction[i] * u;  
212 +    
213 +     fu = model->calcF(tempX);  
214 +
215 +     if(fu <= fMinVar){
216 +
217 +       if(u >= minVar)
218 +         leftVar = minVar;
219 +       else
220 +         rightVar = minVar;
221 +
222 +       v  = prevMinVar;
223 +       prevMinVar = minVar;
224 +       minVar = u;
225 +      
226 +       fv = fPrevMinVar;
227 +       fPrevMinVar = fMinVar;
228 +       fMinVar = fu;
229 +      
230 +     }
231 +     else{
232 +       if (u < minVar) leftVar = u;
233 +         else rightVar= u;
234 +        
235 +       if(fu <= fPrevMinVar || prevMinVar == minVar) {
236 +         v  = prevMinVar;
237 +         fv = fPrevMinVar;
238 +         prevMinVar = u;
239 +         fPrevMinVar = fu;
240 +      }
241 +      else if ( fu <= fv || v == minVar || v == prevMinVar ) {
242 +        v  = u;
243 +         fv = fu;
244 +      }  
245 +    }    
246 +
247 +    midVar = (leftVar + rightVar) /2;
248 +
249 +     if (checkConvergence() > 0){
250 +       minStatus = MINSTATUS_CONVERGE;
251 +       return;
252 +     }
253 +    
254 +  }
255 +
256 +
257    minStatus = MINSTATUS_MAXITER;
258 <  return;
258 >  return;  
259   }
260 +
261 + int BrentMinimizer::checkConvergence(){
262 +  
263 +  if (fabs(minVar - midVar) <  stepTol)
264 +    return 1;
265 +  else
266 +    return -1;
267 + }
268 +
269 + /*******************************************************
270 + *    Bracketing a minimum of a real function Y=F(X)    *
271 + *             using MNBRAK subroutine                  *
272 + * ---------------------------------------------------- *
273 + * REFERENCE: "Numerical recipes, The Art of Scientific *
274 + *             Computing by W.H. Press, B.P. Flannery,  *
275 + *             S.A. Teukolsky and W.T. Vetterling,      *
276 + *             Cambridge university Press, 1986".       *
277 + * ---------------------------------------------------- *
278 + * We have different situation here, we want to limit the
279 + ********************************************************/
280 + void BrentMinimizer::bracket(double& cx, double& fc, double& ax, double& fa, double& bx, double& fb){
281 +  vector<double> currentX;
282 +  vector<double> tempX;
283 +  double u, r, q;
284 +  double fu;
285 +  double ulim;  
286 +  const double TINY = 1.0e-20;
287 +  const double GLIMIT = 100.0;
288 +  const double GoldenRatio = 0.618034;
289 +  const int MAXBRACKETITER = 100;
290 +  currentX = model->getX();
291 +  tempX.resize(currentX.size());
292 +
293 +   if (fb > fa){
294 +     swap(fa, fb);
295 +     swap(ax, bx);
296 +   }
297 +
298 +   cx = bx + GoldenRatio * (bx - ax);
299 +
300 +   fc = model->calcF(tempX);
301 +
302 +   for(int k = 0; k < MAXBRACKETITER && (fb < fc); k++){
303 +    
304 +     r = (bx  - ax) * (fb -fc);
305 +     q = (bx - cx) * (fb - fa);
306 +     u = bx -((bx - cx)*q - (bx-ax)*r)/(2.0 * copysign(max(fabs(q-r), TINY) ,q-r));
307 +     ulim = bx + GLIMIT *(cx - bx);
308 +    
309 +     for (int i = 0; i < tempX.size(); i ++)
310 +       tempX[i] = currentX[i] + direction[i] * u;  
311 +    
312 +     if ((bx -u) * (u -cx) > 0){
313 +       fu = model->calcF(tempX);
314 +
315 +       if (fu < fc){
316 +         ax = bx;
317 +         bx = u;
318 +         fa = fb;
319 +         fb = fu;
320 +       }
321 +       else if (fu > fb){
322 +         cx = u;
323 +         fc = fu;
324 +         return;
325 +       }      
326 +     }
327 +     else if ((cx - u)* (u - ulim) > 0.0){
328 +
329 +       fu = model->calcF(tempX);
330 +
331 +       if (fu < fc){
332 +         bx = cx;
333 +         cx = u;
334 +         u = cx + GoldenRatio * (cx - bx);
335 +
336 +         fb = fc;
337 +         fc = fu;
338 +
339 +         for (int i = 0; i < tempX.size(); i ++)
340 +           tempX[i] = currentX[i] + direction[i] * u;    
341 +
342 +         fu = model->calcF(tempX);
343 +       }
344 +     }
345 +     else if ((u-ulim) * (ulim - cx) >= 0.0){
346 +       u = ulim;
347 +      
348 +       fu =  model->calcF(tempX);
349 +    
350 +     }
351 +     else {
352 +       u = cx + GoldenRatio * (cx -bx);
353 +      
354 +       fu = model->calcF(tempX);
355 +     }
356 +
357 +     ax = bx;
358 +     bx = cx;
359 +     cx = u;
360 +    
361 +     fa = fb;
362 +     fb = fc;
363 +     fc = fu;
364 +
365 +   }
366 +      
367 + }
368 +

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines