OOPSE/libmdtools/Minimizer1D.cpp

#include "Minimizer1D.hpp"
#include "math.h"

//----------------------------------------------------------------------------//
void Minimizer1D::Minimize(vector<double>& direction, double left, double right){
  setDirection(direction);
  setRange(left,right);
  minimize();
}

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

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

  tempX = currentX + direction * alpha;
  fAlpha = model->calcF(tempX);

  tempX = currentX + direction * 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);

      tempX = currentX + beta * direction;

      prevMinVar = minVar;
      fPrevMinVar = fMinVar;
      
      minVar = beta;
      fMinVar = model->calcF(tempX);
      
    }
    else{
      rightVar = beta;
      beta = alpha;
      alpha = leftVar + (1 - goldenRatio) * (rightVar  - leftVar);

      tempX = currentX + alpha * direction;

      prevMinVar = minVar;
      fPrevMinVar = fMinVar;
      
      minVar = alpha;
      fMinVar = model->calcF(tempX);
    }
     
  }

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

  tempX = currentX + leftVar * direction;
  fLeftVar = model->calcF(tempX);
  
  tempX = currentX + rightVar * direction;
  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;
  
  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);

     tempX = currentX + u * direction;
     fu = model->calcF();   

     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;

     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:	1010
Committed:	Tue Feb 3 20:43:08 2004 UTC (20 years, 7 months ago) by tim
File size:	5136 byte(s)
Log Message:	to avoid cyclic depency, refactory constraint class
#	User	Rev	Content
1	tim	1002	#include "Minimizer1D.hpp"
2	tim	1005	#include "math.h"
3
4			//----------------------------------------------------------------------------//
5	tim	1010	void Minimizer1D::Minimize(vector<double>& direction, double left, double right){
6	tim	1002	setDirection(direction);
7			setRange(left,right);
8			minimize();
9			}
10
11	tim	1005	//----------------------------------------------------------------------------//
12			GoldenSectionMinimizer::GoldenSectionMinimizer(NLModel* nlp)
13			:Minimizer1D(nlp){
14			setName("GoldenSection");
15			}
16	tim	1002
17	tim	1005	int GoldenSectionMinimizer::checkConvergence(){
18
19	tim	1002	if ((rightVar - leftVar) < stepTol)
20	tim	1010	return 1;
21	tim	1002	else
22			return -1;
23			}
24			void GoldenSectionMinimizer::minimize(){
25			vector<double> tempX;
26			vector <double> currentX;
27
28			const double goldenRatio = 0.618034;
29
30			currentX = model->getX();
31
32			alpha = leftVar + (1 - goldenRatio) * (rightVar - leftVar);
33			beta = leftVar + goldenRatio * (rightVar - leftVar);
34
35			tempX = currentX + direction * alpha;
36			fAlpha = model->calcF(tempX);
37
38			tempX = currentX + direction * beta;
39			fBeta = model->calcF(tempX);
40
41			for(currentIter = 0; currentIter < maxIteration; currentIter++){
42
43			if (checkConvergence() > 0){
44			minStatus = MINSTATUS_CONVERGE;
45			return;
46			}
47
48			if (fAlpha > fBeta){
49			leftVar = alpha;
50			alpha = beta;
51			beta = leftVar + goldenRatio * (rightVar - leftVar);
52
53			tempX = currentX + beta * direction;
54
55			prevMinVar = minVar;
56			fPrevMinVar = fMinVar;
57
58			minVar = beta;
59			fMinVar = model->calcF(tempX);
60
61			}
62			else{
63			rightVar = beta;
64			beta = alpha;
65			alpha = leftVar + (1 - goldenRatio) * (rightVar - leftVar);
66
67			tempX = currentX + alpha * direction;
68
69			prevMinVar = minVar;
70			fPrevMinVar = fMinVar;
71
72			minVar = alpha;
73			fMinVar = model->calcF(tempX);
74			}
75
76			}
77
78			minStatus = MINSTATUS_MAXITER;
79
80			}
81
82	tim	1005	/**
83			* Brent's method is a root-finding algorithm which combines root bracketing, interval bisection,
84			* and inverse quadratic interpolation.
85	tim	1002	*/
86	tim	1005	BrentMinimizer::BrentMinimizer(NLModel* nlp)
87			:Minimizer1D(nlp){
88			setName("Brent");
89			}
90	tim	1002
91			void BrentMinimizer::minimize(){
92
93	tim	1005	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	tim	1010	double fLeftVar, fRightVar;
101	tim	1005	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	tim	1010	fLeftVar = model->calcF(tempX);
112	tim	1005
113			tempX = currentX + rightVar * direction;
114	tim	1010	fRightVar = model->calcF(tempX);
115	tim	1005
116	tim	1010	if(fRightVar < fLeftVar) {
117			prevMinVar = rightVar;
118			fPrevMinVar = fRightVar;
119	tim	1005	v = leftVar;
120			fv = fLeftVar;
121			}
122			else {
123			prevMinVar = leftVar;
124	tim	1010	fPrevMinVar = fLeftVar;
125	tim	1005	v = rightVar;
126	tim	1010	fv = fRightVar;
127	tim	1005	}
128	tim	1010
129			midVar = leftVar + rightVar;
130	tim	1005
131	tim	1002	for(currentIter = 0; currentIter < maxIteration; currentIter){
132
133	tim	1005	// a trial parabolic fit
134			if (fabs(e) > stepTol){
135	tim	1002
136	tim	1005	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	tim	1002
141	tim	1005	if (q > 0.0)
142			p = -p;
143	tim	1002
144	tim	1005	q = fabs(q);
145
146			etemp = e;
147			e = d;
148
149	tim	1010	if(fabs(p) >= fabs(0.5qetemp) \|\| p <= q(leftVar - minVar) \|\| p >= q(rightVar - minVar)){
150	tim	1005	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	tim	1002	}
209
210
211			minStatus = MINSTATUS_MAXITER;
212	tim	1010	return;
213	tim	1002	}
214	tim	1005
215	tim	1010	int BrentMinimizer::checkConvergence(){
216	tim	1005
217			if (fabs(minVar - midVar) < stepTol)
218			return 1;
219			else
220			return -1;
221			}