OOPSE/libmdtools/Minimizer1D.cpp

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

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;

      prevMinVar = minVar;
      fPrevMinVar = fMinVar;
      
      minVar = beta;
      fMinVar = model->calcF(tempX);
      
    }
    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;

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

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

     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:	1015
Committed:	Tue Feb 3 22:54:52 2004 UTC (20 years, 5 months ago) by tim
File size:	5233 byte(s)
Log Message:	NLModel0, NLModel1 pass uit test
#	Content
1	#include "Minimizer1D.hpp"
2	#include "math.h"
3
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	void GoldenSectionMinimizer::minimize(){
17	vector<double> tempX;
18	vector <double> currentX;
19
20	const double goldenRatio = 0.618034;
21
22	tempX = currentX = model->getX();
23
24	alpha = leftVar + (1 - goldenRatio) * (rightVar - leftVar);
25	beta = leftVar + goldenRatio * (rightVar - leftVar);
26
27	for (int i = 0; i < tempX.size(); i ++)
28	tempX[i] = currentX[i] + direction[i] * alpha;
29
30	fAlpha = model->calcF(tempX);
31
32	for (int i = 0; i < tempX.size(); i ++)
33	tempX[i] = currentX[i] + direction[i] * beta;
34
35	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	beta = leftVar + goldenRatio * (rightVar - leftVar);
48
49	for (int i = 0; i < tempX.size(); i ++)
50	tempX[i] = currentX[i] + direction[i] * beta;
51
52	prevMinVar = minVar;
53	fPrevMinVar = fMinVar;
54
55	minVar = beta;
56	fMinVar = model->calcF(tempX);
57
58	}
59	else{
60	rightVar = beta;
61	beta = alpha;
62	alpha = leftVar + (1 - goldenRatio) * (rightVar - leftVar);
63
64	for (int i = 0; i < tempX.size(); i ++)
65	tempX[i] = currentX[i] + direction[i] * alpha;
66
67	prevMinVar = minVar;
68	fPrevMinVar = fMinVar;
69
70	minVar = alpha;
71	fMinVar = model->calcF(tempX);
72	}
73
74	}
75
76	minStatus = MINSTATUS_MAXITER;
77
78	}
79
80	/**
81	* Brent's method is a root-finding algorithm which combines root bracketing, interval bisection,
82	* and inverse quadratic interpolation.
83	*/
84	BrentMinimizer::BrentMinimizer(NLModel* nlp)
85	:Minimizer1D(nlp){
86	setName("Brent");
87	}
88
89	void BrentMinimizer::minimize(){
90
91	double fu, fv, fw;
92	double p, q, r;
93	double u, v, w;
94	double d;
95	double e;
96	double etemp;
97	double stepTol2;
98	double fLeftVar, fRightVar;
99	const double goldenRatio = 0.3819660;
100	vector<double> tempX, currentX;
101
102	stepTol2 = 2 * stepTol;
103	e = 0;
104	d = 0;
105
106	currentX = tempX = model->getX();
107
108	for (int i = 0; i < tempX.size(); i ++)
109	tempX[i] = currentX[i] + direction[i] * leftVar;
110
111	fLeftVar = model->calcF(tempX);
112
113	for (int i = 0; i < tempX.size(); i ++)
114	tempX[i] = currentX[i] + direction[i] * rightVar;
115
116	fRightVar = model->calcF(tempX);
117
118	if(fRightVar < fLeftVar) {
119	prevMinVar = rightVar;
120	fPrevMinVar = fRightVar;
121	v = leftVar;
122	fv = fLeftVar;
123	}
124	else {
125	prevMinVar = leftVar;
126	fPrevMinVar = fLeftVar;
127	v = rightVar;
128	fv = fRightVar;
129	}
130
131	midVar = leftVar + rightVar;
132
133	for(currentIter = 0; currentIter < maxIteration; currentIter){
134
135	// a trial parabolic fit
136	if (fabs(e) > stepTol){
137
138	r = (minVar - prevMinVar) * (fMinVar - fv);
139	q = (minVar - v) * (fMinVar - fPrevMinVar);
140	p = (minVar - v) q -(minVar - prevMinVar)r;
141	q = 2.0 *(q-r);
142
143	if (q > 0.0)
144	p = -p;
145
146	q = fabs(q);
147
148	etemp = e;
149	e = d;
150
151	if(fabs(p) >= fabs(0.5qetemp) \|\| p <= q(leftVar - minVar) \|\| p >= q(rightVar - minVar)){
152	e = minVar >= midVar ? leftVar - minVar : rightVar - minVar;
153	d = goldenRatio * e;
154	}
155	else{
156	d = p/q;
157	u = minVar + d;
158	if ( u - leftVar < stepTol2 \|\| rightVar - u < stepTol2)
159	d = midVar > minVar ? stepTol : - stepTol;
160	}
161	}
162	//golden section
163	else{
164	e = minVar >=midVar? leftVar - minVar : rightVar - minVar;
165	d =goldenRatio * e;
166	}
167
168	u = fabs(d) >= stepTol ? minVar + d : minVar + copysign(d, stepTol);
169
170	for (int i = 0; i < tempX.size(); i ++)
171	tempX[i] = currentX[i] + direction[i] * u;
172
173	fu = model->calcF();
174
175	if(fu <= fMinVar){
176
177	if(u >= minVar)
178	leftVar = minVar;
179	else
180	rightVar = minVar;
181
182	v = prevMinVar;
183	fv = fPrevMinVar;
184	prevMinVar = minVar;
185	fPrevMinVar = fMinVar;
186	minVar = u;
187	fMinVar = fu;
188
189	}
190	else{
191	if (u < minVar) leftVar = u;
192	else rightVar= u;
193	if(fu <= fPrevMinVar \|\| prevMinVar == minVar) {
194	v = prevMinVar;
195	fv = fPrevMinVar;
196	prevMinVar = u;
197	fPrevMinVar = fu;
198	}
199	else if ( fu <= fv \|\| v == minVar \|\| v == prevMinVar ) {
200	v = u;
201	fv = fu;
202	}
203	}
204
205	midVar = leftVar + rightVar;
206
207	if (checkConvergence() > 0){
208	minStatus = MINSTATUS_CONVERGE;
209	return;
210	}
211
212	}
213
214
215	minStatus = MINSTATUS_MAXITER;
216	return;
217	}
218
219	int BrentMinimizer::checkConvergence(){
220
221	if (fabs(minVar - midVar) < stepTol)
222	return 1;
223	else
224	return -1;
225	}