OOPSE/libmdtools/Minimizer1D.cpp

#include "Minimizer1D.hpp"
#include "math.h"
#include "Utility.hpp"
GoldenSectionMinimizer::GoldenSectionMinimizer(NLModel* nlp)
                               :Minimizer1D(nlp){
  setName("GoldenSection");
}

int GoldenSectionMinimizer::checkConvergence(){

  if ((rightVar - leftVar) < stepTol)
    return 1;
  else
    return -1;
}

void GoldenSectionMinimizer::minimize(){
  vector<double> tempX;
  vector <double> currentX;
  double curF;
  const double goldenRatio = 0.618034;
  
  tempX = currentX =  model->getX();
  model->calcF();
  curF = model->getF();
  
  alpha = leftVar + (1 - goldenRatio) * (rightVar  - leftVar);
  beta = leftVar + goldenRatio * (rightVar - leftVar);

  for (int i = 0; i < tempX.size(); i ++)
    tempX[i] = currentX[i] + direction[i] * alpha;

  fAlpha = model->calcF(tempX);

  for (int i = 0; i < tempX.size(); i ++)
    tempX[i] = currentX[i] + direction[i] * beta;

  fBeta = model->calcF(tempX);

  for(currentIter = 0; currentIter < maxIteration; currentIter++){

     if (checkConvergence() > 0){

      //quick hack
      if (fMinVar > curF) {
        fMinVar = curF;
        minVar = 0;
        minStatus = MINSTATUS_ERROR;
      }
      
       minStatus = MINSTATUS_CONVERGE;
       return;
     }
    
    if (fAlpha > fBeta){
      leftVar = alpha;
      alpha = beta;
      
      beta =  leftVar + goldenRatio * (rightVar - leftVar);

      for (int i = 0; i < tempX.size(); i ++)
        tempX[i] = currentX[i] + direction[i] * beta;
      fAlpha = fBeta;
      fBeta  = model->calcF(tempX);
      
      prevMinVar = alpha;
      fPrevMinVar = fAlpha;      
      minVar = beta;
      fMinVar = fBeta;      
    }
    else{
      rightVar = beta;
      beta = alpha;

      alpha = leftVar + (1 - goldenRatio) * (rightVar  - leftVar);

      for (int i = 0; i < tempX.size(); i ++)
        tempX[i] = currentX[i] + direction[i] * alpha;

      fBeta = fAlpha;
      fAlpha = model->calcF(tempX);

      prevMinVar = beta;
      fPrevMinVar = fBeta;

      minVar = alpha;
      fMinVar = fAlpha;       
    }
     
  }

  cerr << "GoldenSectionMinimizer Warning : can not reach tolerance" << endl;
  minStatus = MINSTATUS_MAXITER;
    
}

/**
 * Brent's method is a root-finding algorithm which combines root bracketing, interval bisection,
 * and inverse quadratic interpolation.
 */
BrentMinimizer::BrentMinimizer(NLModel* nlp)
                   :Minimizer1D(nlp){
  setName("Brent");
}

void BrentMinimizer::minimize(vector<double>& direction, double left, double right){

  //brent algorithm ascending order

  if (left > right)
      setRange(right, left);
  else    
    setRange(left, right);
  
  setDirection(direction);
  
  minimize();
}
void BrentMinimizer::minimize(){

  double fu, fv, fw;
  double p, q, r;
  double u, v, w;
  double d;
  double e;
  double etemp;
  double stepTol2;
  double fLeftVar, fRightVar;
  const double goldenRatio = 0.3819660;
  vector<double> tempX, currentX;
  
  stepTol2 = 2 * stepTol;

  e = 0;
  d = 0;

  currentX = model->getX();
  tempX.resize(currentX.size());

  
  for (int i = 0; i < tempX.size(); i ++)
    tempX[i] = currentX[i] + direction[i] * leftVar;
  
  fLeftVar = model->calcF(tempX);

  for (int i = 0; i < tempX.size(); i ++)
    tempX[i] = currentX[i] + direction[i] * rightVar;  
  
  fRightVar = model->calcF(tempX);

  // find an interior point  left < interior < right which satisfy f(left) > f(interior) and f(right) > f(interior)
   
  bracket(minVar, fMinVar, leftVar, fLeftVar, rightVar, fRightVar);

  if(fRightVar < fLeftVar) {
    prevMinVar = rightVar;
    fPrevMinVar = fRightVar;
    v  = leftVar;
    fv = fLeftVar;
  }
  else {
    prevMinVar = leftVar;
    fPrevMinVar = fLeftVar;
    v  = rightVar;
    fv = fRightVar;
  }
    
  minVar =  rightVar+ goldenRatio * (rightVar - leftVar);

  for (int i = 0; i < tempX.size(); i ++)
    tempX[i] = currentX[i] + direction[i] * minVar; 

  fMinVar = model->calcF(tempX);
  
  prevMinVar = v = minVar;
  fPrevMinVar = fv = fMinVar;
  midVar = (leftVar + rightVar) / 2;
  
  for(currentIter = 0; currentIter < maxIteration; currentIter++){

     //construct a trial parabolic fit
     if (fabs(e) > stepTol){

       r = (minVar - prevMinVar) * (fMinVar - fv);
       q = (minVar - v) * (fMinVar - fPrevMinVar);
       p = (minVar - v) *q -(minVar - prevMinVar)*r;
       q = 2.0 *(q-r);

       if (q > 0.0)
         p = -p;

       q = fabs(q);

       etemp = e;
       e  = d;

       if(fabs(p) >= fabs(0.5*q*etemp) || p <= q*(leftVar - minVar) || p >= q*(rightVar - minVar)){
         //reject parabolic fit and use golden section step instead
         e =  minVar >= midVar ? leftVar - minVar : rightVar - minVar;
         d = goldenRatio * e;
       }
       else{

        //take the parabolic step
         d = p/q;
         u = minVar + d;
         if ( u - leftVar < stepTol2 || rightVar - u  < stepTol2)
           d = midVar > minVar ? stepTol : - stepTol;
       }
       
     }
     else{
       e = minVar >= midVar ? leftVar -minVar : rightVar-minVar;
       d = goldenRatio * e;        
     }

     u = fabs(d) >= stepTol ? minVar + d : minVar + copysign(stepTol, d);

     for (int i = 0; i < tempX.size(); i ++)
       tempX[i] = currentX[i] + direction[i] * u;  
    
     fu = model->calcF(tempX);   

     if(fu <= fMinVar){

       if(u >= minVar)
         leftVar = minVar;
       else
         rightVar = minVar;

       v  = prevMinVar;
       prevMinVar = minVar;
       minVar = u;
       
       fv = fPrevMinVar;
       fPrevMinVar = fMinVar;
       fMinVar = fu;
       
     }
     else{
       if (u < minVar) leftVar = u;
         else rightVar= u;
         
       if(fu <= fPrevMinVar || prevMinVar == minVar) {
         v  = prevMinVar;
         fv = fPrevMinVar;
         prevMinVar = u;
         fPrevMinVar = fu;
      }
      else if ( fu <= fv || v == minVar || v == prevMinVar ) {
        v  = u;
         fv = fu;
      }   
    }    

    midVar = (leftVar + rightVar) /2;

     if (checkConvergence() > 0){
       minStatus = MINSTATUS_CONVERGE;
       return;
     }
    
  }


  minStatus = MINSTATUS_MAXITER;
  return;  
}

int BrentMinimizer::checkConvergence(){
  
  if (fabs(minVar - midVar) <  stepTol)
    return 1;
  else
    return -1;
}

/*******************************************************
*    Bracketing a minimum of a real function Y=F(X)    *
*             using MNBRAK subroutine                  *
* ---------------------------------------------------- *
* REFERENCE: "Numerical recipes, The Art of Scientific *
*             Computing by W.H. Press, B.P. Flannery,  *
*             S.A. Teukolsky and W.T. Vetterling,      *
*             Cambridge university Press, 1986".       *
* ---------------------------------------------------- *
* We have different situation here, we want to limit the 
********************************************************/
void BrentMinimizer::bracket(double& cx, double& fc, double& ax, double& fa, double& bx, double& fb){
  vector<double> currentX;
  vector<double> tempX;
  double u, r, q;
  double fu;
  double ulim;  
  const double TINY = 1.0e-20;
  const double GLIMIT = 100.0;
  const double GoldenRatio = 0.618034;
  const int MAXBRACKETITER = 100; 
  currentX = model->getX();
  tempX.resize(currentX.size());

   if (fb > fa){
     swap(fa, fb);
     swap(ax, bx);
   }

   cx = bx + GoldenRatio * (bx - ax);

   fc = model->calcF(tempX);

   for(int k = 0; k < MAXBRACKETITER && (fb < fc); k++){
    
     r = (bx  - ax) * (fb -fc);
     q = (bx - cx) * (fb - fa);
     u = bx -((bx - cx)*q - (bx-ax)*r)/(2.0 * copysign(max(fabs(q-r), TINY) ,q-r));
     ulim = bx + GLIMIT *(cx - bx);
     
     for (int i = 0; i < tempX.size(); i ++)
       tempX[i] = currentX[i] + direction[i] * u;  
     
     if ((bx -u) * (u -cx) > 0){ 
       fu = model->calcF(tempX);

       if (fu < fc){
         ax = bx;
         bx = u;
         fa = fb;
         fb = fu;
       }
       else if (fu > fb){
         cx = u;
         fc = fu;
         return;
       }       
     }
     else if ((cx - u)* (u - ulim) > 0.0){

       fu = model->calcF(tempX); 

       if (fu < fc){
         bx = cx;
         cx = u;
         u = cx + GoldenRatio * (cx - bx);

         fb = fc;
         fc = fu;

         for (int i = 0; i < tempX.size(); i ++)
           tempX[i] = currentX[i] + direction[i] * u;    

         fu = model->calcF(tempX);
       }
     }
     else if ((u-ulim) * (ulim - cx) >= 0.0){
       u = ulim;
       
       fu =  model->calcF(tempX);
     
     }
     else {
       u = cx + GoldenRatio * (cx -bx);
       
       fu = model->calcF(tempX);
     }

     ax = bx;
     bx = cx;
     cx = u;
     
     fa = fb;
     fb = fc;
     fc = fu;

   }
      
}

Revision:	1064
Committed:	Tue Feb 24 15:44:45 2004 UTC (21 years, 2 months ago) by tim
File size:	9147 byte(s)
Log Message:	Using inherit instead of compose to implement Minimizer both versions are working
#	User	Rev	Content
1	tim	1002	#include "Minimizer1D.hpp"
2	tim	1005	#include "math.h"
3	tim	1023	#include "Utility.hpp"
4	tim	1005	GoldenSectionMinimizer::GoldenSectionMinimizer(NLModel* nlp)
5			:Minimizer1D(nlp){
6			setName("GoldenSection");
7			}
8	tim	1002
9	tim	1005	int GoldenSectionMinimizer::checkConvergence(){
10
11	tim	1002	if ((rightVar - leftVar) < stepTol)
12	tim	1010	return 1;
13	tim	1002	else
14			return -1;
15			}
16	tim	1031
17	tim	1002	void GoldenSectionMinimizer::minimize(){
18			vector<double> tempX;
19			vector <double> currentX;
20	tim	1064	double curF;
21	tim	1002	const double goldenRatio = 0.618034;
22
23	tim	1015	tempX = currentX = model->getX();
24	tim	1064	model->calcF();
25			curF = model->getF();
26
27	tim	1002	alpha = leftVar + (1 - goldenRatio) * (rightVar - leftVar);
28			beta = leftVar + goldenRatio * (rightVar - leftVar);
29
30	tim	1015	for (int i = 0; i < tempX.size(); i ++)
31			tempX[i] = currentX[i] + direction[i] * alpha;
32
33	tim	1002	fAlpha = model->calcF(tempX);
34
35	tim	1015	for (int i = 0; i < tempX.size(); i ++)
36			tempX[i] = currentX[i] + direction[i] * beta;
37
38	tim	1002	fBeta = model->calcF(tempX);
39
40			for(currentIter = 0; currentIter < maxIteration; currentIter++){
41
42			if (checkConvergence() > 0){
43	tim	1064
44			//quick hack
45			if (fMinVar > curF) {
46			fMinVar = curF;
47			minVar = 0;
48			minStatus = MINSTATUS_ERROR;
49			}
50
51	tim	1002	minStatus = MINSTATUS_CONVERGE;
52			return;
53			}
54
55			if (fAlpha > fBeta){
56			leftVar = alpha;
57			alpha = beta;
58	tim	1023
59	tim	1002	beta = leftVar + goldenRatio * (rightVar - leftVar);
60
61	tim	1015	for (int i = 0; i < tempX.size(); i ++)
62			tempX[i] = currentX[i] + direction[i] * beta;
63	tim	1023	fAlpha = fBeta;
64			fBeta = model->calcF(tempX);
65	tim	1002
66	tim	1023	prevMinVar = alpha;
67			fPrevMinVar = fAlpha;
68	tim	1002	minVar = beta;
69	tim	1023	fMinVar = fBeta;
70	tim	1002	}
71			else{
72			rightVar = beta;
73			beta = alpha;
74	tim	1023
75	tim	1002	alpha = leftVar + (1 - goldenRatio) * (rightVar - leftVar);
76
77	tim	1015	for (int i = 0; i < tempX.size(); i ++)
78			tempX[i] = currentX[i] + direction[i] * alpha;
79	tim	1002
80	tim	1023	fBeta = fAlpha;
81			fAlpha = model->calcF(tempX);
82
83			prevMinVar = beta;
84			fPrevMinVar = fBeta;
85
86	tim	1002	minVar = alpha;
87	tim	1023	fMinVar = fAlpha;
88	tim	1002	}
89
90			}
91
92	tim	1057	cerr << "GoldenSectionMinimizer Warning : can not reach tolerance" << endl;
93	tim	1002	minStatus = MINSTATUS_MAXITER;
94
95			}
96
97	tim	1005	/**
98			* Brent's method is a root-finding algorithm which combines root bracketing, interval bisection,
99			* and inverse quadratic interpolation.
100	tim	1002	*/
101	tim	1005	BrentMinimizer::BrentMinimizer(NLModel* nlp)
102			:Minimizer1D(nlp){
103			setName("Brent");
104			}
105	tim	1002
106	tim	1031	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	tim	1002	void BrentMinimizer::minimize(){
120
121	tim	1005	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	tim	1010	double fLeftVar, fRightVar;
129	tim	1005	const double goldenRatio = 0.3819660;
130			vector<double> tempX, currentX;
131
132			stepTol2 = 2 * stepTol;
133	tim	1031
134	tim	1005	e = 0;
135			d = 0;
136
137	tim	1031	currentX = model->getX();
138			tempX.resize(currentX.size());
139	tim	1005
140	tim	1031
141
142
143	tim	1015	for (int i = 0; i < tempX.size(); i ++)
144			tempX[i] = currentX[i] + direction[i] * leftVar;
145
146	tim	1010	fLeftVar = model->calcF(tempX);
147	tim	1015
148			for (int i = 0; i < tempX.size(); i ++)
149			tempX[i] = currentX[i] + direction[i] * rightVar;
150	tim	1005
151	tim	1010	fRightVar = model->calcF(tempX);
152	tim	1005
153	tim	1031	// 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	tim	1010	if(fRightVar < fLeftVar) {
158			prevMinVar = rightVar;
159			fPrevMinVar = fRightVar;
160	tim	1005	v = leftVar;
161			fv = fLeftVar;
162			}
163			else {
164			prevMinVar = leftVar;
165	tim	1010	fPrevMinVar = fLeftVar;
166	tim	1005	v = rightVar;
167	tim	1010	fv = fRightVar;
168	tim	1005	}
169	tim	1031
170			minVar = rightVar+ goldenRatio * (rightVar - leftVar);
171	tim	1010
172	tim	1031	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	tim	1023	midVar = (leftVar + rightVar) / 2;
180	tim	1005
181	tim	1023	for(currentIter = 0; currentIter < maxIteration; currentIter++){
182	tim	1002
183	tim	1031	//construct a trial parabolic fit
184	tim	1005	if (fabs(e) > stepTol){
185	tim	1002
186	tim	1005	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	tim	1002
191	tim	1005	if (q > 0.0)
192			p = -p;
193	tim	1002
194	tim	1005	q = fabs(q);
195
196			etemp = e;
197			e = d;
198
199	tim	1010	if(fabs(p) >= fabs(0.5qetemp) \|\| p <= q(leftVar - minVar) \|\| p >= q(rightVar - minVar)){
200	tim	1031	//reject parabolic fit and use golden section step instead
201	tim	1005	e = minVar >= midVar ? leftVar - minVar : rightVar - minVar;
202			d = goldenRatio * e;
203			}
204			else{
205	tim	1031
206			//take the parabolic step
207	tim	1005	d = p/q;
208			u = minVar + d;
209			if ( u - leftVar < stepTol2 \|\| rightVar - u < stepTol2)
210			d = midVar > minVar ? stepTol : - stepTol;
211			}
212	tim	1031
213	tim	1005	}
214			else{
215	tim	1031	e = minVar >= midVar ? leftVar -minVar : rightVar-minVar;
216			d = goldenRatio * e;
217	tim	1005	}
218
219	tim	1031	u = fabs(d) >= stepTol ? minVar + d : minVar + copysign(stepTol, d);
220	tim	1005
221	tim	1015	for (int i = 0; i < tempX.size(); i ++)
222			tempX[i] = currentX[i] + direction[i] * u;
223
224	tim	1023	fu = model->calcF(tempX);
225	tim	1005
226			if(fu <= fMinVar){
227
228			if(u >= minVar)
229			leftVar = minVar;
230			else
231			rightVar = minVar;
232
233			v = prevMinVar;
234	tim	1031	prevMinVar = minVar;
235			minVar = u;
236
237	tim	1005	fv = fPrevMinVar;
238			fPrevMinVar = fMinVar;
239			fMinVar = fu;
240
241			}
242			else{
243			if (u < minVar) leftVar = u;
244	tim	1031	else rightVar= u;
245
246	tim	1005	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	tim	1023	midVar = (leftVar + rightVar) /2;
259	tim	1005
260			if (checkConvergence() > 0){
261			minStatus = MINSTATUS_CONVERGE;
262			return;
263			}
264
265	tim	1002	}
266
267
268			minStatus = MINSTATUS_MAXITER;
269	tim	1010	return;
270	tim	1002	}
271	tim	1005
272	tim	1010	int BrentMinimizer::checkConvergence(){
273	tim	1005
274			if (fabs(minVar - midVar) < stepTol)
275			return 1;
276			else
277			return -1;
278			}
279	tim	1031
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