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 1031 by tim, Fri Feb 6 18:58:06 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    }
# Line 88 | Line 92 | BrentMinimizer::BrentMinimizer(NLModel* nlp)
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 102 | 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 <  tempX = currentX + leftVar * direction;
129 >  
130 >
131 >  
132 >  for (int i = 0; i < tempX.size(); i ++)
133 >    tempX[i] = currentX[i] + direction[i] * leftVar;
134 >  
135    fLeftVar = model->calcF(tempX);
136 +
137 +  for (int i = 0; i < tempX.size(); i ++)
138 +    tempX[i] = currentX[i] + direction[i] * rightVar;  
139    
113 –  tempX = currentX + rightVar * direction;
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 125 | Line 155 | void BrentMinimizer::minimize(){
155      v  = rightVar;
156      fv = fRightVar;
157    }
158 +    
159 +  minVar =  rightVar+ goldenRatio * (rightVar - leftVar);
160  
161 <  midVar = leftVar + rightVar;
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 <  for(currentIter = 0; currentIter < maxIteration; currentIter){
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 147 | 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       }
160 –     //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 <     tempX = currentX + u * direction;
211 <     fu = model->calcF();  
210 >     for (int i = 0; i < tempX.size(); i ++)
211 >       tempX[i] = currentX[i] + direction[i] * u;  
212 >    
213 >     fu = model->calcF(tempX);  
214  
215       if(fu <= fMinVar){
216  
# Line 176 | Line 220 | void BrentMinimizer::minimize(){
220           rightVar = minVar;
221  
222         v  = prevMinVar;
179 –       fv = fPrevMinVar;
223         prevMinVar = minVar;
181 –       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 198 | Line 244 | void BrentMinimizer::minimize(){
244        }  
245      }    
246  
247 <    midVar = leftVar + rightVar;
247 >    midVar = (leftVar + rightVar) /2;
248  
249       if (checkConvergence() > 0){
250         minStatus = MINSTATUS_CONVERGE;
# Line 219 | Line 265 | int BrentMinimizer::checkConvergence(){
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