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;       
    }
     
  }

  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(){

  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 = tempX = model->getX();

  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);

  if(fRightVar < fLeftVar) {
    prevMinVar = rightVar;
    fPrevMinVar = fRightVar;
    v  = leftVar;
    fv = fLeftVar;
  }
  else {
    prevMinVar = leftVar;
    fPrevMinVar = fLeftVar;
    v  = rightVar;
    fv = fRightVar;
  }

  midVar = (leftVar + rightVar) / 2;
  
  for(currentIter = 0; currentIter < maxIteration; currentIter++){

     // 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)){
         e =  minVar >= midVar ? leftVar - minVar : rightVar - minVar;
         d = goldenRatio * e;
       }
       else{
         d = p/q;
         u = minVar + d;
         if ( u - leftVar < stepTol2 || rightVar - u  < stepTol2)
           d = midVar > minVar ? stepTol : - stepTol;
       }
     }
     //golden section
     else{
       e = minVar >=midVar? leftVar - minVar : rightVar - minVar;
       d =goldenRatio * e;
     }

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

     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;
       fv = fPrevMinVar;
       prevMinVar = minVar;
       fPrevMinVar = fMinVar;
       minVar = u;
       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;
}
Revision:	1023
Committed:	Wed Feb 4 22:26:00 2004 UTC (20 years, 5 months ago) by tim
File size:	5375 byte(s)
Log Message:	Fix a bunch of bugs :-) Single version of conjugate gradient with golden search linesearch pass a couple of functions test. Brent's algorithm is still broken
#	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			void GoldenSectionMinimizer::minimize(){
17			vector<double> tempX;
18			vector <double> currentX;
19
20			const double goldenRatio = 0.618034;
21
22	tim	1015	tempX = currentX = model->getX();
23	tim	1002
24			alpha = leftVar + (1 - goldenRatio) * (rightVar - leftVar);
25			beta = leftVar + goldenRatio * (rightVar - leftVar);
26
27	tim	1015	for (int i = 0; i < tempX.size(); i ++)
28			tempX[i] = currentX[i] + direction[i] * alpha;
29
30	tim	1002	fAlpha = model->calcF(tempX);
31
32	tim	1015	for (int i = 0; i < tempX.size(); i ++)
33			tempX[i] = currentX[i] + direction[i] * beta;
34
35	tim	1002	fBeta = model->calcF(tempX);
36
37			for(currentIter = 0; currentIter < maxIteration; currentIter++){
38
39			if (checkConvergence() > 0){
40			minStatus = MINSTATUS_CONVERGE;
41			return;
42			}
43
44			if (fAlpha > fBeta){
45			leftVar = alpha;
46			alpha = beta;
47	tim	1023
48	tim	1002	beta = leftVar + goldenRatio * (rightVar - leftVar);
49
50	tim	1015	for (int i = 0; i < tempX.size(); i ++)
51			tempX[i] = currentX[i] + direction[i] * beta;
52	tim	1023	fAlpha = fBeta;
53			fBeta = model->calcF(tempX);
54	tim	1002
55	tim	1023	prevMinVar = alpha;
56			fPrevMinVar = fAlpha;
57	tim	1002	minVar = beta;
58	tim	1023	fMinVar = fBeta;
59	tim	1002	}
60			else{
61			rightVar = beta;
62			beta = alpha;
63	tim	1023
64	tim	1002	alpha = leftVar + (1 - goldenRatio) * (rightVar - leftVar);
65
66	tim	1015	for (int i = 0; i < tempX.size(); i ++)
67			tempX[i] = currentX[i] + direction[i] * alpha;
68	tim	1002
69	tim	1023	fBeta = fAlpha;
70			fAlpha = model->calcF(tempX);
71
72			prevMinVar = beta;
73			fPrevMinVar = fBeta;
74
75	tim	1002	minVar = alpha;
76	tim	1023	fMinVar = fAlpha;
77	tim	1002	}
78
79			}
80
81			minStatus = MINSTATUS_MAXITER;
82
83			}
84
85	tim	1005	/**
86			* Brent's method is a root-finding algorithm which combines root bracketing, interval bisection,
87			* and inverse quadratic interpolation.
88	tim	1002	*/
89	tim	1005	BrentMinimizer::BrentMinimizer(NLModel* nlp)
90			:Minimizer1D(nlp){
91			setName("Brent");
92			}
93	tim	1002
94			void BrentMinimizer::minimize(){
95
96	tim	1005	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	tim	1010	double fLeftVar, fRightVar;
104	tim	1005	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	tim	1015	for (int i = 0; i < tempX.size(); i ++)
114			tempX[i] = currentX[i] + direction[i] * leftVar;
115
116	tim	1010	fLeftVar = model->calcF(tempX);
117	tim	1015
118			for (int i = 0; i < tempX.size(); i ++)
119			tempX[i] = currentX[i] + direction[i] * rightVar;
120	tim	1005
121	tim	1010	fRightVar = model->calcF(tempX);
122	tim	1005
123	tim	1010	if(fRightVar < fLeftVar) {
124			prevMinVar = rightVar;
125			fPrevMinVar = fRightVar;
126	tim	1005	v = leftVar;
127			fv = fLeftVar;
128			}
129			else {
130			prevMinVar = leftVar;
131	tim	1010	fPrevMinVar = fLeftVar;
132	tim	1005	v = rightVar;
133	tim	1010	fv = fRightVar;
134	tim	1005	}
135	tim	1010
136	tim	1023	midVar = (leftVar + rightVar) / 2;
137	tim	1005
138	tim	1023	for(currentIter = 0; currentIter < maxIteration; currentIter++){
139	tim	1002
140	tim	1005	// a trial parabolic fit
141			if (fabs(e) > stepTol){
142	tim	1002
143	tim	1005	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	tim	1002
148	tim	1005	if (q > 0.0)
149			p = -p;
150	tim	1002
151	tim	1005	q = fabs(q);
152
153			etemp = e;
154			e = d;
155
156	tim	1010	if(fabs(p) >= fabs(0.5qetemp) \|\| p <= q(leftVar - minVar) \|\| p >= q(rightVar - minVar)){
157	tim	1005	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	tim	1015	for (int i = 0; i < tempX.size(); i ++)
176			tempX[i] = currentX[i] + direction[i] * u;
177
178	tim	1023	fu = model->calcF(tempX);
179	tim	1005
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	tim	1023	midVar = (leftVar + rightVar) /2;
211	tim	1005
212			if (checkConvergence() > 0){
213			minStatus = MINSTATUS_CONVERGE;
214			return;
215			}
216
217	tim	1002	}
218
219
220			minStatus = MINSTATUS_MAXITER;
221	tim	1010	return;
222	tim	1002	}
223	tim	1005
224	tim	1010	int BrentMinimizer::checkConvergence(){
225	tim	1005
226			if (fabs(minVar - midVar) < stepTol)
227			return 1;
228			else
229			return -1;
230			}