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 1023 by tim, Wed Feb 4 22:26:00 2004 UTC vs.
Revision 1064 by tim, Tue Feb 24 15:44:45 2004 UTC

# Line 13 | Line 13 | int GoldenSectionMinimizer::checkConvergence(){
13    else
14      return -1;
15   }
16 +
17   void GoldenSectionMinimizer::minimize(){
18    vector<double> tempX;
19    vector <double> currentX;
20 <
20 >  double curF;
21    const double goldenRatio = 0.618034;
22    
23    tempX = currentX =  model->getX();
24 <
24 >  model->calcF();
25 >  curF = model->getF();
26 >  
27    alpha = leftVar + (1 - goldenRatio) * (rightVar  - leftVar);
28    beta = leftVar + goldenRatio * (rightVar - leftVar);
29  
# Line 37 | Line 40 | void GoldenSectionMinimizer::minimize(){
40    for(currentIter = 0; currentIter < maxIteration; currentIter++){
41  
42       if (checkConvergence() > 0){
43 +
44 +      //quick hack
45 +      if (fMinVar > curF) {
46 +        fMinVar = curF;
47 +        minVar = 0;
48 +        minStatus = MINSTATUS_ERROR;
49 +      }
50 +      
51         minStatus = MINSTATUS_CONVERGE;
52         return;
53       }
# Line 78 | Line 89 | void GoldenSectionMinimizer::minimize(){
89      
90    }
91  
92 +  cerr << "GoldenSectionMinimizer Warning : can not reach tolerance" << endl;
93    minStatus = MINSTATUS_MAXITER;
94      
95   }
# Line 91 | Line 103 | BrentMinimizer::BrentMinimizer(NLModel* nlp)
103    setName("Brent");
104   }
105  
106 + void BrentMinimizer::minimize(vector<double>& direction, double left, double right){
107 +
108 +  //brent algorithm ascending order
109 +
110 +  if (left > right)
111 +      setRange(right, left);
112 +  else    
113 +    setRange(left, right);
114 +  
115 +  setDirection(direction);
116 +  
117 +  minimize();
118 + }
119   void BrentMinimizer::minimize(){
120  
121    double fu, fv, fw;
# Line 105 | Line 130 | void BrentMinimizer::minimize(){
130    vector<double> tempX, currentX;
131    
132    stepTol2 = 2 * stepTol;
133 +
134    e = 0;
135    d = 0;
136  
137 <  currentX = tempX = model->getX();
137 >  currentX = model->getX();
138 >  tempX.resize(currentX.size());
139  
140 +  
141 +
142 +  
143    for (int i = 0; i < tempX.size(); i ++)
144      tempX[i] = currentX[i] + direction[i] * leftVar;
145    
# Line 120 | Line 150 | void BrentMinimizer::minimize(){
150    
151    fRightVar = model->calcF(tempX);
152  
153 +  // find an interior point  left < interior < right which satisfy f(left) > f(interior) and f(right) > f(interior)
154 +  
155 +  bracket(minVar, fMinVar, leftVar, fLeftVar, rightVar, fRightVar);
156 +
157    if(fRightVar < fLeftVar) {
158      prevMinVar = rightVar;
159      fPrevMinVar = fRightVar;
# Line 132 | Line 166 | void BrentMinimizer::minimize(){
166      v  = rightVar;
167      fv = fRightVar;
168    }
169 +    
170 +  minVar =  rightVar+ goldenRatio * (rightVar - leftVar);
171  
172 +  for (int i = 0; i < tempX.size(); i ++)
173 +    tempX[i] = currentX[i] + direction[i] * minVar;
174 +
175 +  fMinVar = model->calcF(tempX);
176 +  
177 +  prevMinVar = v = minVar;
178 +  fPrevMinVar = fv = fMinVar;
179    midVar = (leftVar + rightVar) / 2;
180    
181    for(currentIter = 0; currentIter < maxIteration; currentIter++){
182  
183 <     // a trial parabolic fit
183 >     //construct a trial parabolic fit
184       if (fabs(e) > stepTol){
185  
186         r = (minVar - prevMinVar) * (fMinVar - fv);
# Line 154 | Line 197 | void BrentMinimizer::minimize(){
197         e  = d;
198  
199         if(fabs(p) >= fabs(0.5*q*etemp) || p <= q*(leftVar - minVar) || p >= q*(rightVar - minVar)){
200 +         //reject parabolic fit and use golden section step instead
201           e =  minVar >= midVar ? leftVar - minVar : rightVar - minVar;
202           d = goldenRatio * e;
203         }
204         else{
205 +
206 +        //take the parabolic step
207           d = p/q;
208           u = minVar + d;
209           if ( u - leftVar < stepTol2 || rightVar - u  < stepTol2)
210             d = midVar > minVar ? stepTol : - stepTol;
211         }
212 +      
213       }
167 –     //golden section
214       else{
215 <       e = minVar >=midVar? leftVar - minVar : rightVar - minVar;
216 <       d =goldenRatio * e;
215 >       e = minVar >= midVar ? leftVar -minVar : rightVar-minVar;
216 >       d = goldenRatio * e;        
217       }
218  
219 <     u = fabs(d) >= stepTol ? minVar + d : minVar + copysign(d, stepTol);
219 >     u = fabs(d) >= stepTol ? minVar + d : minVar + copysign(stepTol, d);
220  
221       for (int i = 0; i < tempX.size(); i ++)
222         tempX[i] = currentX[i] + direction[i] * u;  
# Line 185 | Line 231 | void BrentMinimizer::minimize(){
231           rightVar = minVar;
232  
233         v  = prevMinVar;
188 –       fv = fPrevMinVar;
234         prevMinVar = minVar;
190 –       fPrevMinVar = fMinVar;
235         minVar = u;
236 +      
237 +       fv = fPrevMinVar;
238 +       fPrevMinVar = fMinVar;
239         fMinVar = fu;
240        
241       }
242       else{
243         if (u < minVar) leftVar = u;
244 <       else rightVar= u;
244 >         else rightVar= u;
245 >        
246         if(fu <= fPrevMinVar || prevMinVar == minVar) {
247           v  = prevMinVar;
248           fv = fPrevMinVar;
# Line 228 | Line 276 | int BrentMinimizer::checkConvergence(){
276    else
277      return -1;
278   }
279 +
280 + /*******************************************************
281 + *    Bracketing a minimum of a real function Y=F(X)    *
282 + *             using MNBRAK subroutine                  *
283 + * ---------------------------------------------------- *
284 + * REFERENCE: "Numerical recipes, The Art of Scientific *
285 + *             Computing by W.H. Press, B.P. Flannery,  *
286 + *             S.A. Teukolsky and W.T. Vetterling,      *
287 + *             Cambridge university Press, 1986".       *
288 + * ---------------------------------------------------- *
289 + * We have different situation here, we want to limit the
290 + ********************************************************/
291 + void BrentMinimizer::bracket(double& cx, double& fc, double& ax, double& fa, double& bx, double& fb){
292 +  vector<double> currentX;
293 +  vector<double> tempX;
294 +  double u, r, q;
295 +  double fu;
296 +  double ulim;  
297 +  const double TINY = 1.0e-20;
298 +  const double GLIMIT = 100.0;
299 +  const double GoldenRatio = 0.618034;
300 +  const int MAXBRACKETITER = 100;
301 +  currentX = model->getX();
302 +  tempX.resize(currentX.size());
303 +
304 +   if (fb > fa){
305 +     swap(fa, fb);
306 +     swap(ax, bx);
307 +   }
308 +
309 +   cx = bx + GoldenRatio * (bx - ax);
310 +
311 +   fc = model->calcF(tempX);
312 +
313 +   for(int k = 0; k < MAXBRACKETITER && (fb < fc); k++){
314 +    
315 +     r = (bx  - ax) * (fb -fc);
316 +     q = (bx - cx) * (fb - fa);
317 +     u = bx -((bx - cx)*q - (bx-ax)*r)/(2.0 * copysign(max(fabs(q-r), TINY) ,q-r));
318 +     ulim = bx + GLIMIT *(cx - bx);
319 +    
320 +     for (int i = 0; i < tempX.size(); i ++)
321 +       tempX[i] = currentX[i] + direction[i] * u;  
322 +    
323 +     if ((bx -u) * (u -cx) > 0){
324 +       fu = model->calcF(tempX);
325 +
326 +       if (fu < fc){
327 +         ax = bx;
328 +         bx = u;
329 +         fa = fb;
330 +         fb = fu;
331 +       }
332 +       else if (fu > fb){
333 +         cx = u;
334 +         fc = fu;
335 +         return;
336 +       }      
337 +     }
338 +     else if ((cx - u)* (u - ulim) > 0.0){
339 +
340 +       fu = model->calcF(tempX);
341 +
342 +       if (fu < fc){
343 +         bx = cx;
344 +         cx = u;
345 +         u = cx + GoldenRatio * (cx - bx);
346 +
347 +         fb = fc;
348 +         fc = fu;
349 +
350 +         for (int i = 0; i < tempX.size(); i ++)
351 +           tempX[i] = currentX[i] + direction[i] * u;    
352 +
353 +         fu = model->calcF(tempX);
354 +       }
355 +     }
356 +     else if ((u-ulim) * (ulim - cx) >= 0.0){
357 +       u = ulim;
358 +      
359 +       fu =  model->calcF(tempX);
360 +    
361 +     }
362 +     else {
363 +       u = cx + GoldenRatio * (cx -bx);
364 +      
365 +       fu = model->calcF(tempX);
366 +     }
367 +
368 +     ax = bx;
369 +     bx = cx;
370 +     cx = u;
371 +    
372 +     fa = fb;
373 +     fb = fc;
374 +     fc = fu;
375 +
376 +   }
377 +      
378 + }
379 +

Diff Legend

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