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 fLeft, fRight;
  const double goldenRatio = 0.3819660;
  vector<double> tempX, currentX;
  
  stepTol2 = 2 * stepTol;
  e = 0;
  d = 0;

  currentX = tempX = model->getX();

  tempX = currentX + leftVar * direction;
  fLeft = model->calcF(tempX);
  
  tempX = currentX + rightVar * direction;
  fRight = model->calcF(tempX);

  if(fRight < fLeft) {
    prevMinPoint = rightVar;
    fPrevMinVar = fRight;
    v  = leftVar;
    fv = fLeftVar;
  }
  else {
    prevMinVar = leftVar;
    fPrevMinVar = fLeft;
    v  = rightVar;
    fv = fRight;
  }
  
  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;

//----------------------------------------------------------------------------//
  
}

BrentMinimizer::checkConvergence(){
  
  if (fabs(minVar - midVar) <  stepTol)
    return 1;
  else
    return -1;
}
Revision:	1005
Committed:	Tue Feb 3 15:21:32 2004 UTC (20 years, 5 months ago) by tim
File size:	5162 byte(s)
Log Message:	begin unit test of minimizer
#	User	Rev	Content
1	tim	1002	#include "Minimizer1D.hpp"
2	tim	1005	#include "math.h"
3
4			//----------------------------------------------------------------------------//
5	tim	1002	void Minimizer1D::Minimize(vector<double>& direction), double left, double right); {
6			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	1005	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			double fLeft, fRight;
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			fLeft = model->calcF(tempX);
112
113			tempX = currentX + rightVar * direction;
114			fRight = model->calcF(tempX);
115
116			if(fRight < fLeft) {
117			prevMinPoint = rightVar;
118			fPrevMinVar = fRight;
119			v = leftVar;
120			fv = fLeftVar;
121			}
122			else {
123			prevMinVar = leftVar;
124			fPrevMinVar = fLeft;
125			v = rightVar;
126			fv = fRight;
127			}
128
129	tim	1002	for(currentIter = 0; currentIter < maxIteration; currentIter){
130
131	tim	1005	// a trial parabolic fit
132			if (fabs(e) > stepTol){
133	tim	1002
134	tim	1005	r = (minVar - prevMinVar) * (fMinVar - fv);
135			q = (minVar - v) * (fMinVar - fPrevMinVar);
136			p = (minVar - v) q -(minVar - prevMinVar)r;
137			q = 2.0 *(q-r);
138	tim	1002
139	tim	1005	if (q > 0.0)
140			p = -p;
141	tim	1002
142	tim	1005	q = fabs(q);
143
144			etemp = e;
145			e = d;
146
147			if(fabs(p) >= fabs(0.5qetemp)) \|\| p <= q(leftVar - minVar) \|\| p >= q(rightVar - minVar)){
148			e = minVar >= midVar ? leftVar - minVar : rightVar - minVar;
149			d = goldenRatio * e;
150			}
151			else{
152			d = p/q;
153			u = minVar + d;
154			if ( u - leftVar < stepTol2 \|\| rightVar - u < stepTol2)
155			d = midVar > minVar ? stepTol : - stepTol;
156			}
157			}
158			//golden section
159			else{
160			e = minVar >=midVar? leftVar - minVar : rightVar - minVar;
161			d =goldenRatio * e;
162			}
163
164			u = fabs(d) >= stepTol ? minVar + d : minVar + copysign(d, stepTol);
165
166			tempX = currentX + u * direction;
167			fu = model->calcF();
168
169			if(fu <= fMinVar){
170
171			if(u >= minVar)
172			leftVar = minVar;
173			else
174			rightVar = minVar;
175
176			v = prevMinVar;
177			fv = fPrevMinVar;
178			prevMinVar = minVar;
179			fPrevMinVar = fMinVar;
180			minVar = u;
181			fMinVar = fu;
182
183			}
184			else{
185			if (u < minVar) leftVar = u;
186			else rightVar= u;
187			if(fu <= fPrevMinVar \|\| prevMinVar == minVar) {
188			v = prevMinVar;
189			fv = fPrevMinVar;
190			prevMinVar = u;
191			fPrevMinVar = fu;
192			}
193			else if ( fu <= fv \|\| v == minVar \|\| v == prevMinVar ) {
194			v = u;
195			fv = fu;
196			}
197			}
198
199			midVar = leftVar + rightVar;
200
201			if (checkConvergence() > 0){
202			minStatus = MINSTATUS_CONVERGE;
203			return;
204			}
205
206	tim	1002	}
207
208
209			minStatus = MINSTATUS_MAXITER;
210			return;
211	tim	1005
212			//----------------------------------------------------------------------------//
213
214	tim	1002	}
215	tim	1005
216			BrentMinimizer::checkConvergence(){
217
218			if (fabs(minVar - midVar) < stepTol)
219			return 1;
220			else
221			return -1;
222			}