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 1031 by tim, Fri Feb 6 18:58:06 2004 UTC

# Line 13 | Line 13 | void GoldenSectionMinimizer::minimize(){
13    else
14      return -1;
15   }
16 +
17   void GoldenSectionMinimizer::minimize(){
18    vector<double> tempX;
19    vector <double> currentX;
# Line 91 | Line 92 | void BrentMinimizer::minimize(){
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    double fu, fv, fw;
# Line 105 | Line 119 | void BrentMinimizer::minimize(){
119    vector<double> tempX, currentX;
120    
121    stepTol2 = 2 * stepTol;
122 +
123    e = 0;
124    d = 0;
125  
126 <  currentX = tempX = model->getX();
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    
# Line 120 | Line 139 | void BrentMinimizer::minimize(){
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;
# Line 132 | Line 155 | void BrentMinimizer::minimize(){
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 <     // a trial parabolic fit
172 >     //construct a trial parabolic fit
173       if (fabs(e) > stepTol){
174  
175         r = (minVar - prevMinVar) * (fMinVar - fv);
# Line 154 | Line 186 | void BrentMinimizer::minimize(){
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       }
167 –     //golden section
203       else{
204 <       e = minVar >=midVar? leftVar - minVar : rightVar - minVar;
205 <       d =goldenRatio * e;
204 >       e = minVar >= midVar ? leftVar -minVar : rightVar-minVar;
205 >       d = goldenRatio * e;        
206       }
207  
208 <     u = fabs(d) >= stepTol ? minVar + d : minVar + copysign(d, stepTol);
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;  
# Line 185 | Line 220 | void BrentMinimizer::minimize(){
220           rightVar = minVar;
221  
222         v  = prevMinVar;
188 –       fv = fPrevMinVar;
223         prevMinVar = minVar;
190 –       fPrevMinVar = fMinVar;
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;
233 >         else rightVar= u;
234 >        
235         if(fu <= fPrevMinVar || prevMinVar == minVar) {
236           v  = prevMinVar;
237           fv = fPrevMinVar;
# Line 227 | Line 264 | int BrentMinimizer::checkConvergence(){
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