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 1002 by tim, Mon Feb 2 20:29:41 2004 UTC vs.
Revision 1064 by tim, Tue Feb 24 15:44:45 2004 UTC

# Line 1 | Line 1
1   #include "Minimizer1D.hpp"
2 < void Minimizer1D::Minimize(vector<double>& direction), double left, double right); {
3 <  setDirection(direction);
4 <  setRange(left,right);
5 <  minimize();
2 > #include "math.h"
3 > #include "Utility.hpp"
4 > GoldenSectionMinimizer::GoldenSectionMinimizer(NLModel* nlp)
5 >                               :Minimizer1D(nlp){
6 >  setName("GoldenSection");
7   }
8  
9 < int Minimizer1D::checkConvergence(){
9 > int GoldenSectionMinimizer::checkConvergence(){
10  
11    if ((rightVar - leftVar) < stepTol)
12 <    return
12 >    return 1;
13    else
14      return -1;
15   }
# Line 16 | Line 17 | void GoldenSectionMinimizer::minimize(){
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 40 | 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);
52 <      
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   }
96  
97 < /*
98 < *
97 > /**
98 > * Brent's method is a root-finding algorithm which combines root bracketing, interval bisection,
99 > * and inverse quadratic interpolation.
100   */
101 + BrentMinimizer::BrentMinimizer(NLModel* nlp)
102 +                   :Minimizer1D(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 <  for(currentIter = 0; currentIter < maxIteration; currentIter){
121 >  double fu, fv, fw;
122 >  double p, q, r;
123 >  double u, v, w;
124 >  double d;
125 >  double e;
126 >  double etemp;
127 >  double stepTol2;
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 = model->getX();
138 +  tempX.resize(currentX.size());
139  
140 +  
141  
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 = fLeftVar;
166 +    v  = rightVar;
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 +  prevMinVar = v = minVar;
178 +  fPrevMinVar = fv = fMinVar;
179 +  midVar = (leftVar + rightVar) / 2;
180 +  
181 +  for(currentIter = 0; currentIter < maxIteration; currentIter++){
182 +
183 +     //construct a trial parabolic fit
184 +     if (fabs(e) > stepTol){
185 +
186 +       r = (minVar - prevMinVar) * (fMinVar - fv);
187 +       q = (minVar - v) * (fMinVar - fPrevMinVar);
188 +       p = (minVar - v) *q -(minVar - prevMinVar)*r;
189 +       q = 2.0 *(q-r);
190 +
191 +       if (q > 0.0)
192 +         p = -p;
193 +
194 +       q = fabs(q);
195 +
196 +       etemp = e;
197 +       e  = d;
198 +
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 +     }
214 +     else{
215 +       e = minVar >= midVar ? leftVar -minVar : rightVar-minVar;
216 +       d = goldenRatio * e;        
217 +     }
218 +
219 +     u = fabs(d) >= stepTol ? minVar + d : minVar + copysign(stepTol, d);
220 +
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 +
228 +       if(u >= minVar)
229 +         leftVar = minVar;
230 +       else
231 +         rightVar = minVar;
232 +
233 +       v  = prevMinVar;
234 +       prevMinVar = minVar;
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;
245 +        
246 +       if(fu <= fPrevMinVar || prevMinVar == minVar) {
247 +         v  = prevMinVar;
248 +         fv = fPrevMinVar;
249 +         prevMinVar = u;
250 +         fPrevMinVar = fu;
251 +      }
252 +      else if ( fu <= fv || v == minVar || v == prevMinVar ) {
253 +        v  = u;
254 +         fv = fu;
255 +      }  
256 +    }    
257 +
258 +    midVar = (leftVar + rightVar) /2;
259 +
260 +     if (checkConvergence() > 0){
261 +       minStatus = MINSTATUS_CONVERGE;
262 +       return;
263 +     }
264 +    
265 +  }
266 +
267 +
268    minStatus = MINSTATUS_MAXITER;
269 <  return;
269 >  return;  
270   }
271 +
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