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 1023 by tim, Wed Feb 4 22:26:00 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   }
15 –
16   void GoldenSectionMinimizer::minimize(){
17    vector<double> tempX;
18    vector <double> currentX;
19  
20    const double goldenRatio = 0.618034;
21    
22 <  currentX =  model->getX();
22 >  tempX = currentX =  model->getX();
23  
24    alpha = leftVar + (1 - goldenRatio) * (rightVar  - leftVar);
25    beta = leftVar + goldenRatio * (rightVar - leftVar);
26  
27 <  tempX = currentX + direction * alpha;
27 >  for (int i = 0; i < tempX.size(); i ++)
28 >    tempX[i] = currentX[i] + direction[i] * alpha;
29 >
30    fAlpha = model->calcF(tempX);
31  
32 <  tempX = currentX + direction * beta;
32 >  for (int i = 0; i < tempX.size(); i ++)
33 >    tempX[i] = currentX[i] + direction[i] * beta;
34 >
35    fBeta = model->calcF(tempX);
36  
37    for(currentIter = 0; currentIter < maxIteration; currentIter++){
# Line 40 | Line 44 | void GoldenSectionMinimizer::minimize(){
44      if (fAlpha > fBeta){
45        leftVar = alpha;
46        alpha = beta;
47 +      
48        beta =  leftVar + goldenRatio * (rightVar - leftVar);
49  
50 <      tempX = currentX + beta * direction;
51 <
52 <      prevMinVar = minVar;
53 <      fPrevMinVar = fMinVar;
50 >      for (int i = 0; i < tempX.size(); i ++)
51 >        tempX[i] = currentX[i] + direction[i] * beta;
52 >      fAlpha = fBeta;
53 >      fBeta  = model->calcF(tempX);
54        
55 +      prevMinVar = alpha;
56 +      fPrevMinVar = fAlpha;      
57        minVar = beta;
58 <      fMinVar = model->calcF(tempX);
52 <      
58 >      fMinVar = fBeta;      
59      }
60      else{
61        rightVar = beta;
62        beta = alpha;
63 +
64        alpha = leftVar + (1 - goldenRatio) * (rightVar  - leftVar);
65  
66 <      tempX = currentX + alpha * direction;
66 >      for (int i = 0; i < tempX.size(); i ++)
67 >        tempX[i] = currentX[i] + direction[i] * alpha;
68  
69 <      prevMinVar = minVar;
70 <      fPrevMinVar = fMinVar;
71 <      
69 >      fBeta = fAlpha;
70 >      fAlpha = model->calcF(tempX);
71 >
72 >      prevMinVar = beta;
73 >      fPrevMinVar = fBeta;
74 >
75        minVar = alpha;
76 <      fMinVar = model->calcF(tempX);
76 >      fMinVar = fAlpha;      
77      }
78      
79    }
# Line 71 | Line 82 | void GoldenSectionMinimizer::minimize(){
82      
83   }
84  
85 < /*
86 < *
85 > /**
86 > * Brent's method is a root-finding algorithm which combines root bracketing, interval bisection,
87 > * and inverse quadratic interpolation.
88   */
89 + BrentMinimizer::BrentMinimizer(NLModel* nlp)
90 +                   :Minimizer1D(nlp){
91 +  setName("Brent");
92 + }
93  
94   void BrentMinimizer::minimize(){
95  
96 <  for(currentIter = 0; currentIter < maxIteration; currentIter){
96 >  double fu, fv, fw;
97 >  double p, q, r;
98 >  double u, v, w;
99 >  double d;
100 >  double e;
101 >  double etemp;
102 >  double stepTol2;
103 >  double fLeftVar, fRightVar;
104 >  const double goldenRatio = 0.3819660;
105 >  vector<double> tempX, currentX;
106 >  
107 >  stepTol2 = 2 * stepTol;
108 >  e = 0;
109 >  d = 0;
110  
111 +  currentX = tempX = model->getX();
112  
113 +  for (int i = 0; i < tempX.size(); i ++)
114 +    tempX[i] = currentX[i] + direction[i] * leftVar;
115 +  
116 +  fLeftVar = model->calcF(tempX);
117  
118 +  for (int i = 0; i < tempX.size(); i ++)
119 +    tempX[i] = currentX[i] + direction[i] * rightVar;  
120 +  
121 +  fRightVar = model->calcF(tempX);
122  
123 +  if(fRightVar < fLeftVar) {
124 +    prevMinVar = rightVar;
125 +    fPrevMinVar = fRightVar;
126 +    v  = leftVar;
127 +    fv = fLeftVar;
128    }
129 +  else {
130 +    prevMinVar = leftVar;
131 +    fPrevMinVar = fLeftVar;
132 +    v  = rightVar;
133 +    fv = fRightVar;
134 +  }
135  
136 +  midVar = (leftVar + rightVar) / 2;
137 +  
138 +  for(currentIter = 0; currentIter < maxIteration; currentIter++){
139  
140 +     // a trial parabolic fit
141 +     if (fabs(e) > stepTol){
142 +
143 +       r = (minVar - prevMinVar) * (fMinVar - fv);
144 +       q = (minVar - v) * (fMinVar - fPrevMinVar);
145 +       p = (minVar - v) *q -(minVar - prevMinVar)*r;
146 +       q = 2.0 *(q-r);
147 +
148 +       if (q > 0.0)
149 +         p = -p;
150 +
151 +       q = fabs(q);
152 +
153 +       etemp = e;
154 +       e  = d;
155 +
156 +       if(fabs(p) >= fabs(0.5*q*etemp) || p <= q*(leftVar - minVar) || p >= q*(rightVar - minVar)){
157 +         e =  minVar >= midVar ? leftVar - minVar : rightVar - minVar;
158 +         d = goldenRatio * e;
159 +       }
160 +       else{
161 +         d = p/q;
162 +         u = minVar + d;
163 +         if ( u - leftVar < stepTol2 || rightVar - u  < stepTol2)
164 +           d = midVar > minVar ? stepTol : - stepTol;
165 +       }
166 +     }
167 +     //golden section
168 +     else{
169 +       e = minVar >=midVar? leftVar - minVar : rightVar - minVar;
170 +       d =goldenRatio * e;
171 +     }
172 +
173 +     u = fabs(d) >= stepTol ? minVar + d : minVar + copysign(d, stepTol);
174 +
175 +     for (int i = 0; i < tempX.size(); i ++)
176 +       tempX[i] = currentX[i] + direction[i] * u;  
177 +    
178 +     fu = model->calcF(tempX);  
179 +
180 +     if(fu <= fMinVar){
181 +
182 +       if(u >= minVar)
183 +         leftVar = minVar;
184 +       else
185 +         rightVar = minVar;
186 +
187 +       v  = prevMinVar;
188 +       fv = fPrevMinVar;
189 +       prevMinVar = minVar;
190 +       fPrevMinVar = fMinVar;
191 +       minVar = u;
192 +       fMinVar = fu;
193 +      
194 +     }
195 +     else{
196 +       if (u < minVar) leftVar = u;
197 +       else rightVar= u;
198 +       if(fu <= fPrevMinVar || prevMinVar == minVar) {
199 +         v  = prevMinVar;
200 +         fv = fPrevMinVar;
201 +         prevMinVar = u;
202 +         fPrevMinVar = fu;
203 +      }
204 +      else if ( fu <= fv || v == minVar || v == prevMinVar ) {
205 +        v  = u;
206 +         fv = fu;
207 +      }  
208 +    }    
209 +
210 +    midVar = (leftVar + rightVar) /2;
211 +
212 +     if (checkConvergence() > 0){
213 +       minStatus = MINSTATUS_CONVERGE;
214 +       return;
215 +     }
216 +    
217 +  }
218 +
219 +
220    minStatus = MINSTATUS_MAXITER;
221 <  return;
221 >  return;  
222   }
223 +
224 + int BrentMinimizer::checkConvergence(){
225 +  
226 +  if (fabs(minVar - midVar) <  stepTol)
227 +    return 1;
228 +  else
229 +    return -1;
230 + }

Diff Legend

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