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 1005 by tim, Tue Feb 3 15:21:32 2004 UTC vs.
Revision 1064 by tim, Tue Feb 24 15:44:45 2004 UTC

# Line 1 | Line 1
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 < //----------------------------------------------------------------------------//
3 > #include "Utility.hpp"
4   GoldenSectionMinimizer::GoldenSectionMinimizer(NLModel* nlp)
5                                 :Minimizer1D(nlp){
6    setName("GoldenSection");
# Line 17 | Line 9 | int GoldenSectionMinimizer::checkConvergence(){
9   int GoldenSectionMinimizer::checkConvergence(){
10  
11    if ((rightVar - leftVar) < stepTol)
12 <    return 1
12 >    return 1;
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 <  currentX =  model->getX();
24 <
23 >  tempX = currentX =  model->getX();
24 >  model->calcF();
25 >  curF = model->getF();
26 >  
27    alpha = leftVar + (1 - goldenRatio) * (rightVar  - leftVar);
28    beta = leftVar + goldenRatio * (rightVar - leftVar);
29  
30 <  tempX = currentX + direction * alpha;
30 >  for (int i = 0; i < tempX.size(); i ++)
31 >    tempX[i] = currentX[i] + direction[i] * alpha;
32 >
33    fAlpha = model->calcF(tempX);
34  
35 <  tempX = currentX + direction * beta;
35 >  for (int i = 0; i < tempX.size(); i ++)
36 >    tempX[i] = currentX[i] + direction[i] * beta;
37 >
38    fBeta = model->calcF(tempX);
39  
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 48 | Line 55 | void GoldenSectionMinimizer::minimize(){
55      if (fAlpha > fBeta){
56        leftVar = alpha;
57        alpha = beta;
58 +      
59        beta =  leftVar + goldenRatio * (rightVar - leftVar);
60  
61 <      tempX = currentX + beta * direction;
62 <
63 <      prevMinVar = minVar;
64 <      fPrevMinVar = fMinVar;
61 >      for (int i = 0; i < tempX.size(); i ++)
62 >        tempX[i] = currentX[i] + direction[i] * beta;
63 >      fAlpha = fBeta;
64 >      fBeta  = model->calcF(tempX);
65        
66 +      prevMinVar = alpha;
67 +      fPrevMinVar = fAlpha;      
68        minVar = beta;
69 <      fMinVar = model->calcF(tempX);
60 <      
69 >      fMinVar = fBeta;      
70      }
71      else{
72        rightVar = beta;
73        beta = alpha;
74 +
75        alpha = leftVar + (1 - goldenRatio) * (rightVar  - leftVar);
76  
77 <      tempX = currentX + alpha * direction;
77 >      for (int i = 0; i < tempX.size(); i ++)
78 >        tempX[i] = currentX[i] + direction[i] * alpha;
79  
80 <      prevMinVar = minVar;
81 <      fPrevMinVar = fMinVar;
82 <      
80 >      fBeta = fAlpha;
81 >      fAlpha = model->calcF(tempX);
82 >
83 >      prevMinVar = beta;
84 >      fPrevMinVar = fBeta;
85 >
86        minVar = alpha;
87 <      fMinVar = model->calcF(tempX);
87 >      fMinVar = fAlpha;      
88      }
89      
90    }
91  
92 +  cerr << "GoldenSectionMinimizer Warning : can not reach tolerance" << endl;
93    minStatus = MINSTATUS_MAXITER;
94      
95   }
# Line 88 | 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 97 | Line 125 | void BrentMinimizer::minimize(){
125    double e;
126    double etemp;
127    double stepTol2;
128 <  double fLeft, fRight;
128 >  double fLeftVar, fRightVar;
129    const double goldenRatio = 0.3819660;
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  
110 –  tempX = currentX + leftVar * direction;
111 –  fLeft = model->calcF(tempX);
140    
113 –  tempX = currentX + rightVar * direction;
114 –  fRight = model->calcF(tempX);
141  
142 <  if(fRight < fLeft) {
143 <    prevMinPoint = rightVar;
144 <    fPrevMinVar = fRight;
142 >  
143 >  for (int i = 0; i < tempX.size(); i ++)
144 >    tempX[i] = currentX[i] + direction[i] * leftVar;
145 >  
146 >  fLeftVar = model->calcF(tempX);
147 >
148 >  for (int i = 0; i < tempX.size(); i ++)
149 >    tempX[i] = currentX[i] + direction[i] * rightVar;  
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;
160      v  = leftVar;
161      fv = fLeftVar;
162    }
163    else {
164      prevMinVar = leftVar;
165 <    fPrevMinVar = fLeft;
165 >    fPrevMinVar = fLeftVar;
166      v  = rightVar;
167 <    fv = fRight;
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 <  for(currentIter = 0; currentIter < maxIteration; currentIter){
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 144 | Line 196 | void BrentMinimizer::minimize(){
196         etemp = e;
197         e  = d;
198  
199 <       if(fabs(p) >= fabs(0.5*q*etemp)) || p <= q*(leftVar - minVar) || p >= q*(rightVar - minVar)){
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       }
158 –     //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 <     tempX = currentX + u * direction;
222 <     fu = model->calcF();  
221 >     for (int i = 0; i < tempX.size(); i ++)
222 >       tempX[i] = currentX[i] + direction[i] * u;  
223 >    
224 >     fu = model->calcF(tempX);  
225  
226       if(fu <= fMinVar){
227  
# Line 174 | Line 231 | void BrentMinimizer::minimize(){
231           rightVar = minVar;
232  
233         v  = prevMinVar;
177 –       fv = fPrevMinVar;
234         prevMinVar = minVar;
179 –       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 196 | Line 255 | void BrentMinimizer::minimize(){
255        }  
256      }    
257  
258 <    midVar = leftVar + rightVar;
258 >    midVar = (leftVar + rightVar) /2;
259  
260       if (checkConvergence() > 0){
261         minStatus = MINSTATUS_CONVERGE;
# Line 207 | Line 266 | void BrentMinimizer::minimize(){
266  
267  
268    minStatus = MINSTATUS_MAXITER;
269 <  return;
211 <
212 < //----------------------------------------------------------------------------//
213 <  
269 >  return;  
270   }
271  
272 < BrentMinimizer::checkConvergence(){
272 > int BrentMinimizer::checkConvergence(){
273    
274    if (fabs(minVar - midVar) <  stepTol)
275      return 1;
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