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 1057 by tim, Tue Feb 17 19:23:44 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 78 | Line 79 | void GoldenSectionMinimizer::minimize(){
79      
80    }
81  
82 +  cerr << "GoldenSectionMinimizer Warning : can not reach tolerance" << endl;
83    minStatus = MINSTATUS_MAXITER;
84      
85   }
# Line 91 | Line 93 | void BrentMinimizer::minimize(){
93    setName("Brent");
94   }
95  
96 + void BrentMinimizer::minimize(vector<double>& direction, double left, double right){
97 +
98 +  //brent algorithm ascending order
99 +
100 +  if (left > right)
101 +      setRange(right, left);
102 +  else    
103 +    setRange(left, right);
104 +  
105 +  setDirection(direction);
106 +  
107 +  minimize();
108 + }
109   void BrentMinimizer::minimize(){
110  
111    double fu, fv, fw;
# Line 105 | Line 120 | void BrentMinimizer::minimize(){
120    vector<double> tempX, currentX;
121    
122    stepTol2 = 2 * stepTol;
123 +
124    e = 0;
125    d = 0;
126  
127 <  currentX = tempX = model->getX();
127 >  currentX = model->getX();
128 >  tempX.resize(currentX.size());
129  
130 +  
131 +
132 +  
133    for (int i = 0; i < tempX.size(); i ++)
134      tempX[i] = currentX[i] + direction[i] * leftVar;
135    
# Line 120 | Line 140 | void BrentMinimizer::minimize(){
140    
141    fRightVar = model->calcF(tempX);
142  
143 +  // find an interior point  left < interior < right which satisfy f(left) > f(interior) and f(right) > f(interior)
144 +  
145 +  bracket(minVar, fMinVar, leftVar, fLeftVar, rightVar, fRightVar);
146 +
147    if(fRightVar < fLeftVar) {
148      prevMinVar = rightVar;
149      fPrevMinVar = fRightVar;
# Line 132 | Line 156 | void BrentMinimizer::minimize(){
156      v  = rightVar;
157      fv = fRightVar;
158    }
159 +    
160 +  minVar =  rightVar+ goldenRatio * (rightVar - leftVar);
161  
162 +  for (int i = 0; i < tempX.size(); i ++)
163 +    tempX[i] = currentX[i] + direction[i] * minVar;
164 +
165 +  fMinVar = model->calcF(tempX);
166 +  
167 +  prevMinVar = v = minVar;
168 +  fPrevMinVar = fv = fMinVar;
169    midVar = (leftVar + rightVar) / 2;
170    
171    for(currentIter = 0; currentIter < maxIteration; currentIter++){
172  
173 <     // a trial parabolic fit
173 >     //construct a trial parabolic fit
174       if (fabs(e) > stepTol){
175  
176         r = (minVar - prevMinVar) * (fMinVar - fv);
# Line 154 | Line 187 | void BrentMinimizer::minimize(){
187         e  = d;
188  
189         if(fabs(p) >= fabs(0.5*q*etemp) || p <= q*(leftVar - minVar) || p >= q*(rightVar - minVar)){
190 +         //reject parabolic fit and use golden section step instead
191           e =  minVar >= midVar ? leftVar - minVar : rightVar - minVar;
192           d = goldenRatio * e;
193         }
194         else{
195 +
196 +        //take the parabolic step
197           d = p/q;
198           u = minVar + d;
199           if ( u - leftVar < stepTol2 || rightVar - u  < stepTol2)
200             d = midVar > minVar ? stepTol : - stepTol;
201         }
202 +      
203       }
167 –     //golden section
204       else{
205 <       e = minVar >=midVar? leftVar - minVar : rightVar - minVar;
206 <       d =goldenRatio * e;
205 >       e = minVar >= midVar ? leftVar -minVar : rightVar-minVar;
206 >       d = goldenRatio * e;        
207       }
208  
209 <     u = fabs(d) >= stepTol ? minVar + d : minVar + copysign(d, stepTol);
209 >     u = fabs(d) >= stepTol ? minVar + d : minVar + copysign(stepTol, d);
210  
211       for (int i = 0; i < tempX.size(); i ++)
212         tempX[i] = currentX[i] + direction[i] * u;  
# Line 185 | Line 221 | void BrentMinimizer::minimize(){
221           rightVar = minVar;
222  
223         v  = prevMinVar;
188 –       fv = fPrevMinVar;
224         prevMinVar = minVar;
190 –       fPrevMinVar = fMinVar;
225         minVar = u;
226 +      
227 +       fv = fPrevMinVar;
228 +       fPrevMinVar = fMinVar;
229         fMinVar = fu;
230        
231       }
232       else{
233         if (u < minVar) leftVar = u;
234 <       else rightVar= u;
234 >         else rightVar= u;
235 >        
236         if(fu <= fPrevMinVar || prevMinVar == minVar) {
237           v  = prevMinVar;
238           fv = fPrevMinVar;
# Line 228 | Line 266 | int BrentMinimizer::checkConvergence(){
266    else
267      return -1;
268   }
269 +
270 + /*******************************************************
271 + *    Bracketing a minimum of a real function Y=F(X)    *
272 + *             using MNBRAK subroutine                  *
273 + * ---------------------------------------------------- *
274 + * REFERENCE: "Numerical recipes, The Art of Scientific *
275 + *             Computing by W.H. Press, B.P. Flannery,  *
276 + *             S.A. Teukolsky and W.T. Vetterling,      *
277 + *             Cambridge university Press, 1986".       *
278 + * ---------------------------------------------------- *
279 + * We have different situation here, we want to limit the
280 + ********************************************************/
281 + void BrentMinimizer::bracket(double& cx, double& fc, double& ax, double& fa, double& bx, double& fb){
282 +  vector<double> currentX;
283 +  vector<double> tempX;
284 +  double u, r, q;
285 +  double fu;
286 +  double ulim;  
287 +  const double TINY = 1.0e-20;
288 +  const double GLIMIT = 100.0;
289 +  const double GoldenRatio = 0.618034;
290 +  const int MAXBRACKETITER = 100;
291 +  currentX = model->getX();
292 +  tempX.resize(currentX.size());
293 +
294 +   if (fb > fa){
295 +     swap(fa, fb);
296 +     swap(ax, bx);
297 +   }
298 +
299 +   cx = bx + GoldenRatio * (bx - ax);
300 +
301 +   fc = model->calcF(tempX);
302 +
303 +   for(int k = 0; k < MAXBRACKETITER && (fb < fc); k++){
304 +    
305 +     r = (bx  - ax) * (fb -fc);
306 +     q = (bx - cx) * (fb - fa);
307 +     u = bx -((bx - cx)*q - (bx-ax)*r)/(2.0 * copysign(max(fabs(q-r), TINY) ,q-r));
308 +     ulim = bx + GLIMIT *(cx - bx);
309 +    
310 +     for (int i = 0; i < tempX.size(); i ++)
311 +       tempX[i] = currentX[i] + direction[i] * u;  
312 +    
313 +     if ((bx -u) * (u -cx) > 0){
314 +       fu = model->calcF(tempX);
315 +
316 +       if (fu < fc){
317 +         ax = bx;
318 +         bx = u;
319 +         fa = fb;
320 +         fb = fu;
321 +       }
322 +       else if (fu > fb){
323 +         cx = u;
324 +         fc = fu;
325 +         return;
326 +       }      
327 +     }
328 +     else if ((cx - u)* (u - ulim) > 0.0){
329 +
330 +       fu = model->calcF(tempX);
331 +
332 +       if (fu < fc){
333 +         bx = cx;
334 +         cx = u;
335 +         u = cx + GoldenRatio * (cx - bx);
336 +
337 +         fb = fc;
338 +         fc = fu;
339 +
340 +         for (int i = 0; i < tempX.size(); i ++)
341 +           tempX[i] = currentX[i] + direction[i] * u;    
342 +
343 +         fu = model->calcF(tempX);
344 +       }
345 +     }
346 +     else if ((u-ulim) * (ulim - cx) >= 0.0){
347 +       u = ulim;
348 +      
349 +       fu =  model->calcF(tempX);
350 +    
351 +     }
352 +     else {
353 +       u = cx + GoldenRatio * (cx -bx);
354 +      
355 +       fu = model->calcF(tempX);
356 +     }
357 +
358 +     ax = bx;
359 +     bx = cx;
360 +     cx = u;
361 +    
362 +     fa = fb;
363 +     fb = fc;
364 +     fc = fu;
365 +
366 +   }
367 +      
368 + }
369 +

Diff Legend

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