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 1010 by tim, Tue Feb 3 20:43:08 2004 UTC

# Line 1 | Line 1
1   #include "Minimizer1D.hpp"
2 < void Minimizer1D::Minimize(vector<double>& direction), double left, double right); {
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 < int Minimizer1D::checkConvergence(){
11 > //----------------------------------------------------------------------------//
12 > GoldenSectionMinimizer::GoldenSectionMinimizer(NLModel* nlp)
13 >                               :Minimizer1D(nlp){
14 >  setName("GoldenSection");
15 > }
16  
17 + int GoldenSectionMinimizer::checkConvergence(){
18 +
19    if ((rightVar - leftVar) < stepTol)
20 <    return
20 >    return 1;
21    else
22      return -1;
23   }
15
24   void GoldenSectionMinimizer::minimize(){
25    vector<double> tempX;
26    vector <double> currentX;
# Line 71 | Line 79 | void GoldenSectionMinimizer::minimize(){
79      
80   }
81  
82 < /*
83 < *
82 > /**
83 > * Brent's method is a root-finding algorithm which combines root bracketing, interval bisection,
84 > * and inverse quadratic interpolation.
85   */
86 + BrentMinimizer::BrentMinimizer(NLModel* nlp)
87 +                   :Minimizer1D(nlp){
88 +  setName("Brent");
89 + }
90  
91   void BrentMinimizer::minimize(){
92  
93 +  double fu, fv, fw;
94 +  double p, q, r;
95 +  double u, v, w;
96 +  double d;
97 +  double e;
98 +  double etemp;
99 +  double stepTol2;
100 +  double fLeftVar, fRightVar;
101 +  const double goldenRatio = 0.3819660;
102 +  vector<double> tempX, currentX;
103 +  
104 +  stepTol2 = 2 * stepTol;
105 +  e = 0;
106 +  d = 0;
107 +
108 +  currentX = tempX = model->getX();
109 +
110 +  tempX = currentX + leftVar * direction;
111 +  fLeftVar = model->calcF(tempX);
112 +  
113 +  tempX = currentX + rightVar * direction;
114 +  fRightVar = model->calcF(tempX);
115 +
116 +  if(fRightVar < fLeftVar) {
117 +    prevMinVar = rightVar;
118 +    fPrevMinVar = fRightVar;
119 +    v  = leftVar;
120 +    fv = fLeftVar;
121 +  }
122 +  else {
123 +    prevMinVar = leftVar;
124 +    fPrevMinVar = fLeftVar;
125 +    v  = rightVar;
126 +    fv = fRightVar;
127 +  }
128 +
129 +  midVar = leftVar + rightVar;
130 +  
131    for(currentIter = 0; currentIter < maxIteration; currentIter){
132  
133 +     // a trial parabolic fit
134 +     if (fabs(e) > stepTol){
135  
136 +       r = (minVar - prevMinVar) * (fMinVar - fv);
137 +       q = (minVar - v) * (fMinVar - fPrevMinVar);
138 +       p = (minVar - v) *q -(minVar - prevMinVar)*r;
139 +       q = 2.0 *(q-r);
140  
141 +       if (q > 0.0)
142 +         p = -p;
143  
144 +       q = fabs(q);
145 +
146 +       etemp = e;
147 +       e  = d;
148 +
149 +       if(fabs(p) >= fabs(0.5*q*etemp) || p <= q*(leftVar - minVar) || p >= q*(rightVar - minVar)){
150 +         e =  minVar >= midVar ? leftVar - minVar : rightVar - minVar;
151 +         d = goldenRatio * e;
152 +       }
153 +       else{
154 +         d = p/q;
155 +         u = minVar + d;
156 +         if ( u - leftVar < stepTol2 || rightVar - u  < stepTol2)
157 +           d = midVar > minVar ? stepTol : - stepTol;
158 +       }
159 +     }
160 +     //golden section
161 +     else{
162 +       e = minVar >=midVar? leftVar - minVar : rightVar - minVar;
163 +       d =goldenRatio * e;
164 +     }
165 +
166 +     u = fabs(d) >= stepTol ? minVar + d : minVar + copysign(d, stepTol);
167 +
168 +     tempX = currentX + u * direction;
169 +     fu = model->calcF();  
170 +
171 +     if(fu <= fMinVar){
172 +
173 +       if(u >= minVar)
174 +         leftVar = minVar;
175 +       else
176 +         rightVar = minVar;
177 +
178 +       v  = prevMinVar;
179 +       fv = fPrevMinVar;
180 +       prevMinVar = minVar;
181 +       fPrevMinVar = fMinVar;
182 +       minVar = u;
183 +       fMinVar = fu;
184 +      
185 +     }
186 +     else{
187 +       if (u < minVar) leftVar = u;
188 +       else rightVar= u;
189 +       if(fu <= fPrevMinVar || prevMinVar == minVar) {
190 +         v  = prevMinVar;
191 +         fv = fPrevMinVar;
192 +         prevMinVar = u;
193 +         fPrevMinVar = fu;
194 +      }
195 +      else if ( fu <= fv || v == minVar || v == prevMinVar ) {
196 +        v  = u;
197 +         fv = fu;
198 +      }  
199 +    }    
200 +
201 +    midVar = leftVar + rightVar;
202 +
203 +     if (checkConvergence() > 0){
204 +       minStatus = MINSTATUS_CONVERGE;
205 +       return;
206 +     }
207 +    
208    }
209  
210  
211    minStatus = MINSTATUS_MAXITER;
212 <  return;
212 >  return;  
213   }
214 +
215 + int BrentMinimizer::checkConvergence(){
216 +  
217 +  if (fabs(minVar - midVar) <  stepTol)
218 +    return 1;
219 +  else
220 +    return -1;
221 + }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines