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;

  const double goldenRatio = 0.618034;
  
  tempX = currentX =  model->getX();

  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){
       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:	1057
Committed:	Tue Feb 17 19:23:44 2004 UTC (21 years, 2 months ago) by tim
File size:	8935 byte(s)
Log Message:	adding function shakeF in order to remove the constraint force along bond direction
#	Content
1	#include "Minimizer1D.hpp"
2	#include "math.h"
3	#include "Utility.hpp"
4	GoldenSectionMinimizer::GoldenSectionMinimizer(NLModel* nlp)
5	:Minimizer1D(nlp){
6	setName("GoldenSection");
7	}
8
9	int GoldenSectionMinimizer::checkConvergence(){
10
11	if ((rightVar - leftVar) < stepTol)
12	return 1;
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	tempX = currentX = model->getX();
24
25	alpha = leftVar + (1 - goldenRatio) * (rightVar - leftVar);
26	beta = leftVar + goldenRatio * (rightVar - leftVar);
27
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	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++){
39
40	if (checkConvergence() > 0){
41	minStatus = MINSTATUS_CONVERGE;
42	return;
43	}
44
45	if (fAlpha > fBeta){
46	leftVar = alpha;
47	alpha = beta;
48
49	beta = leftVar + goldenRatio * (rightVar - leftVar);
50
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 = fBeta;
60	}
61	else{
62	rightVar = beta;
63	beta = alpha;
64
65	alpha = leftVar + (1 - goldenRatio) * (rightVar - leftVar);
66
67	for (int i = 0; i < tempX.size(); i ++)
68	tempX[i] = currentX[i] + direction[i] * alpha;
69
70	fBeta = fAlpha;
71	fAlpha = model->calcF(tempX);
72
73	prevMinVar = beta;
74	fPrevMinVar = fBeta;
75
76	minVar = alpha;
77	fMinVar = fAlpha;
78	}
79
80	}
81
82	cerr << "GoldenSectionMinimizer Warning : can not reach tolerance" << endl;
83	minStatus = MINSTATUS_MAXITER;
84
85	}
86
87	/**
88	* Brent's method is a root-finding algorithm which combines root bracketing, interval bisection,
89	* and inverse quadratic interpolation.
90	*/
91	BrentMinimizer::BrentMinimizer(NLModel* nlp)
92	:Minimizer1D(nlp){
93	setName("Brent");
94	}
95
96	void BrentMinimizer::minimize(vector<double>& direction, double left, double right){
97
98	//brent algorithm ascending order
99
100	if (left > right)
101	setRange(right, left);
102	else
103	setRange(left, right);
104
105	setDirection(direction);
106
107	minimize();
108	}
109	void BrentMinimizer::minimize(){
110
111	double fu, fv, fw;
112	double p, q, r;
113	double u, v, w;
114	double d;
115	double e;
116	double etemp;
117	double stepTol2;
118	double fLeftVar, fRightVar;
119	const double goldenRatio = 0.3819660;
120	vector<double> tempX, currentX;
121
122	stepTol2 = 2 * stepTol;
123
124	e = 0;
125	d = 0;
126
127	currentX = model->getX();
128	tempX.resize(currentX.size());
129
130
131
132
133	for (int i = 0; i < tempX.size(); i ++)
134	tempX[i] = currentX[i] + direction[i] * leftVar;
135
136	fLeftVar = model->calcF(tempX);
137
138	for (int i = 0; i < tempX.size(); i ++)
139	tempX[i] = currentX[i] + direction[i] * rightVar;
140
141	fRightVar = model->calcF(tempX);
142
143	// find an interior point left < interior < right which satisfy f(left) > f(interior) and f(right) > f(interior)
144
145	bracket(minVar, fMinVar, leftVar, fLeftVar, rightVar, fRightVar);
146
147	if(fRightVar < fLeftVar) {
148	prevMinVar = rightVar;
149	fPrevMinVar = fRightVar;
150	v = leftVar;
151	fv = fLeftVar;
152	}
153	else {
154	prevMinVar = leftVar;
155	fPrevMinVar = fLeftVar;
156	v = rightVar;
157	fv = fRightVar;
158	}
159
160	minVar = rightVar+ goldenRatio * (rightVar - leftVar);
161
162	for (int i = 0; i < tempX.size(); i ++)
163	tempX[i] = currentX[i] + direction[i] * minVar;
164
165	fMinVar = model->calcF(tempX);
166
167	prevMinVar = v = minVar;
168	fPrevMinVar = fv = fMinVar;
169	midVar = (leftVar + rightVar) / 2;
170
171	for(currentIter = 0; currentIter < maxIteration; currentIter++){
172
173	//construct a trial parabolic fit
174	if (fabs(e) > stepTol){
175
176	r = (minVar - prevMinVar) * (fMinVar - fv);
177	q = (minVar - v) * (fMinVar - fPrevMinVar);
178	p = (minVar - v) q -(minVar - prevMinVar)r;
179	q = 2.0 *(q-r);
180
181	if (q > 0.0)
182	p = -p;
183
184	q = fabs(q);
185
186	etemp = e;
187	e = d;
188
189	if(fabs(p) >= fabs(0.5qetemp) \|\| p <= q(leftVar - minVar) \|\| p >= q(rightVar - minVar)){
190	//reject parabolic fit and use golden section step instead
191	e = minVar >= midVar ? leftVar - minVar : rightVar - minVar;
192	d = goldenRatio * e;
193	}
194	else{
195
196	//take the parabolic step
197	d = p/q;
198	u = minVar + d;
199	if ( u - leftVar < stepTol2 \|\| rightVar - u < stepTol2)
200	d = midVar > minVar ? stepTol : - stepTol;
201	}
202
203	}
204	else{
205	e = minVar >= midVar ? leftVar -minVar : rightVar-minVar;
206	d = goldenRatio * e;
207	}
208
209	u = fabs(d) >= stepTol ? minVar + d : minVar + copysign(stepTol, d);
210
211	for (int i = 0; i < tempX.size(); i ++)
212	tempX[i] = currentX[i] + direction[i] * u;
213
214	fu = model->calcF(tempX);
215
216	if(fu <= fMinVar){
217
218	if(u >= minVar)
219	leftVar = minVar;
220	else
221	rightVar = minVar;
222
223	v = prevMinVar;
224	prevMinVar = minVar;
225	minVar = u;
226
227	fv = fPrevMinVar;
228	fPrevMinVar = fMinVar;
229	fMinVar = fu;
230
231	}
232	else{
233	if (u < minVar) leftVar = u;
234	else rightVar= u;
235
236	if(fu <= fPrevMinVar \|\| prevMinVar == minVar) {
237	v = prevMinVar;
238	fv = fPrevMinVar;
239	prevMinVar = u;
240	fPrevMinVar = fu;
241	}
242	else if ( fu <= fv \|\| v == minVar \|\| v == prevMinVar ) {
243	v = u;
244	fv = fu;
245	}
246	}
247
248	midVar = (leftVar + rightVar) /2;
249
250	if (checkConvergence() > 0){
251	minStatus = MINSTATUS_CONVERGE;
252	return;
253	}
254
255	}
256
257
258	minStatus = MINSTATUS_MAXITER;
259	return;
260	}
261
262	int BrentMinimizer::checkConvergence(){
263
264	if (fabs(minVar - midVar) < stepTol)
265	return 1;
266	else
267	return -1;
268	}
269
270	/*******************************************************
271	* Bracketing a minimum of a real function Y=F(X) *
272	* using MNBRAK subroutine *
273	* ---------------------------------------------------- *
274	* REFERENCE: "Numerical recipes, The Art of Scientific *
275	* Computing by W.H. Press, B.P. Flannery, *
276	* S.A. Teukolsky and W.T. Vetterling, *
277	* Cambridge university Press, 1986". *
278	* ---------------------------------------------------- *
279	* We have different situation here, we want to limit the
280	********************************************************/
281	void BrentMinimizer::bracket(double& cx, double& fc, double& ax, double& fa, double& bx, double& fb){
282	vector<double> currentX;
283	vector<double> tempX;
284	double u, r, q;
285	double fu;
286	double ulim;
287	const double TINY = 1.0e-20;
288	const double GLIMIT = 100.0;
289	const double GoldenRatio = 0.618034;
290	const int MAXBRACKETITER = 100;
291	currentX = model->getX();
292	tempX.resize(currentX.size());
293
294	if (fb > fa){
295	swap(fa, fb);
296	swap(ax, bx);
297	}
298
299	cx = bx + GoldenRatio * (bx - ax);
300
301	fc = model->calcF(tempX);
302
303	for(int k = 0; k < MAXBRACKETITER && (fb < fc); k++){
304
305	r = (bx - ax) * (fb -fc);
306	q = (bx - cx) * (fb - fa);
307	u = bx -((bx - cx)q - (bx-ax)r)/(2.0 * copysign(max(fabs(q-r), TINY) ,q-r));
308	ulim = bx + GLIMIT *(cx - bx);
309
310	for (int i = 0; i < tempX.size(); i ++)
311	tempX[i] = currentX[i] + direction[i] * u;
312
313	if ((bx -u) * (u -cx) > 0){
314	fu = model->calcF(tempX);
315
316	if (fu < fc){
317	ax = bx;
318	bx = u;
319	fa = fb;
320	fb = fu;
321	}
322	else if (fu > fb){
323	cx = u;
324	fc = fu;
325	return;
326	}
327	}
328	else if ((cx - u)* (u - ulim) > 0.0){
329
330	fu = model->calcF(tempX);
331
332	if (fu < fc){
333	bx = cx;
334	cx = u;
335	u = cx + GoldenRatio * (cx - bx);
336
337	fb = fc;
338	fc = fu;
339
340	for (int i = 0; i < tempX.size(); i ++)
341	tempX[i] = currentX[i] + direction[i] * u;
342
343	fu = model->calcF(tempX);
344	}
345	}
346	else if ((u-ulim) * (ulim - cx) >= 0.0){
347	u = ulim;
348
349	fu = model->calcF(tempX);
350
351	}
352	else {
353	u = cx + GoldenRatio * (cx -bx);
354
355	fu = model->calcF(tempX);
356	}
357
358	ax = bx;
359	bx = cx;
360	cx = u;
361
362	fa = fb;
363	fb = fc;
364	fc = fu;
365
366	}
367
368	}
369