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 1010 by tim, Tue Feb 3 20:43:08 2004 UTC vs.
Revision 1057 by tim, Tue Feb 17 19:23:44 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 21 | 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  
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 48 | 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);
60 <      
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    }
81  
82 +  cerr << "GoldenSectionMinimizer Warning : can not reach tolerance" << endl;
83    minStatus = MINSTATUS_MAXITER;
84      
85   }
# Line 88 | Line 93 | BrentMinimizer::BrentMinimizer(NLModel* nlp)
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 102 | 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 <  tempX = currentX + leftVar * direction;
130 >  
131 >
132 >  
133 >  for (int i = 0; i < tempX.size(); i ++)
134 >    tempX[i] = currentX[i] + direction[i] * leftVar;
135 >  
136    fLeftVar = model->calcF(tempX);
137 +
138 +  for (int i = 0; i < tempX.size(); i ++)
139 +    tempX[i] = currentX[i] + direction[i] * rightVar;  
140    
113 –  tempX = currentX + rightVar * direction;
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 125 | Line 156 | void BrentMinimizer::minimize(){
156      v  = rightVar;
157      fv = fRightVar;
158    }
159 +    
160 +  minVar =  rightVar+ goldenRatio * (rightVar - leftVar);
161  
162 <  midVar = leftVar + rightVar;
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 <  for(currentIter = 0; currentIter < maxIteration; currentIter){
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 147 | 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       }
160 –     //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 <     tempX = currentX + u * direction;
212 <     fu = model->calcF();  
211 >     for (int i = 0; i < tempX.size(); i ++)
212 >       tempX[i] = currentX[i] + direction[i] * u;  
213 >    
214 >     fu = model->calcF(tempX);  
215  
216       if(fu <= fMinVar){
217  
# Line 176 | Line 221 | void BrentMinimizer::minimize(){
221           rightVar = minVar;
222  
223         v  = prevMinVar;
179 –       fv = fPrevMinVar;
224         prevMinVar = minVar;
181 –       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 198 | Line 245 | void BrentMinimizer::minimize(){
245        }  
246      }    
247  
248 <    midVar = leftVar + rightVar;
248 >    midVar = (leftVar + rightVar) /2;
249  
250       if (checkConvergence() > 0){
251         minStatus = MINSTATUS_CONVERGE;
# Line 219 | 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