OOPSE/libmdtools/ConjugateMinimizer.cpp

#include "ConjugateMinimizer.hpp"
#include "Utility.hpp"

bool ConjugateMinimizerBase::isSolvable(){

  //conjuage gradient can only solve unconstrained nonlinear model

  if (!model->hasConstraints())
    return true;
  else
    return false;
}

void ConjugateMinimizerBase::Init(){

}

void ConjugateMinimizerBase::printMinizerInfo(){

}

void ConjugateMinimizerBase::Minimize(){
  int maxIteration;
  int nextResetIter;
  int resetFrq;
  int nextWriteIter;
  int writeFrq;
  int lsStatus;
  double gamma;
  double lamda;
  
  
  if (!isSolvable()){
    cout << "ConjugateMinimizerBase Error: This nonlinear model can not be solved by " << minimizerName <<endl;

    exit(1);
  }

  printMinizerInfo();

  resetFrq = paramSet->getResetFrq();
  nextResetIter = resetFrq;

  writeFrq = paramSet->getWriteFrq();
  nextWriteIter = writeFrq;
  
  prevGrad = model->calcGrad();

  direction = prevGrad;
  
  maxIteration = paramSet->getMaxIteration();

  for(currentIter = 0;currentIter < maxIteration; currentIter++){

    // perform line search to minimize f(x + stepSize * direction) where stepSize > 0
    lsMinimizer->minimize(direction, 0.0, 1.0);
    
    lsStatus = lsMinimizer->getMinimizationStatus();
    
    if(lsStatus ==MINSTATUS_ERROR){
      minStatus = MINSTATUS_ERROR;
      return;
    }
    
    prevMinX = minX;
    lamda = lsMinimizer->getMinVar();

    for(int i = 0; i < direction.size(); i++)
      minX[i] = minX[i] + lamda * direction[i];

    //calculate the gradient
    prevGrad = gradient;
    
    gradient = model->calcGrad();

    // stop if converge
    if (checkConvergence() > 0){
      writeOut(minX, currentIter);

      minStatus = MINSTATUS_CONVERGE;
      return;
    }
          

    //calculate the 
    gamma = calcGamma(gradient, prevGrad);

    // update new direction
    prevDirection = direction;

    for(int i = 0; i < direction.size(); i++)  
      direction[i] += gamma * direction[i];

    //
    if (currentIter == nextWriteIter){
      nextWriteIter += writeFrq;
      writeOut(minX, currentIter); 
    }

    if (currentIter == nextResetIter){
      reset();
      nextResetIter += resetFrq;      
    }

  }

  // if writeFrq is not a multipiler of maxIteration, we need to write the final result
  // otherwise, we already write it inside the loop, just skip it
  if(currentIter != (nextWriteIter - writeFrq))
    writeOut(minX, currentIter);

  minStatus = MINSTATUS_MAXITER;
  return;
}

int ConjugateMinimizerBase::checkConvergence(){

  //test absolute gradient tolerance
  
  if (sqrt(dot(gradient, gradient)) < paramSet->getGradTol())
    return 1;
  else
    return -1;
}

void ConjugateMinimizerBase::reset(){  
  
}

double FRCGMinimizer::calcGamma(vector<double>& newGrad, vector<double>& oldGrad){
  return dot(newGrad, newGrad) / dot(oldGrad, newGrad);
}

double PRCGMinimizer::calcGamma(vector<double>& newGrad, vector<double>& oldGrad){
  double gamma;
  vector<double> deltaGrad;
  
  for(int i = 0; i < newGrad.size(); i++)
    deltaGrad.push_back(newGrad[i] - oldGrad[i]);

  return dot(deltaGrad, newGrad) / dot(oldGrad, oldGrad);
  
}
Revision:	1015
Committed:	Tue Feb 3 22:54:52 2004 UTC (20 years, 7 months ago) by tim
File size:	3255 byte(s)
Log Message:	NLModel0, NLModel1 pass uit test
#	User	Rev	Content
1	tim	1011	#include "ConjugateMinimizer.hpp"
2	tim	1015	#include "Utility.hpp"
3
4	tim	1011	bool ConjugateMinimizerBase::isSolvable(){
5
6			//conjuage gradient can only solve unconstrained nonlinear model
7
8	tim	1015	if (!model->hasConstraints())
9	tim	1011	return true;
10			else
11			return false;
12			}
13
14			void ConjugateMinimizerBase::Init(){
15
16			}
17
18			void ConjugateMinimizerBase::printMinizerInfo(){
19
20			}
21
22			void ConjugateMinimizerBase::Minimize(){
23			int maxIteration;
24	tim	1015	int nextResetIter;
25	tim	1011	int resetFrq;
26			int nextWriteIter;
27			int writeFrq;
28	tim	1015	int lsStatus;
29			double gamma;
30			double lamda;
31	tim	1011
32	tim	1015
33	tim	1011	if (!isSolvable()){
34	tim	1015	cout << "ConjugateMinimizerBase Error: This nonlinear model can not be solved by " << minimizerName <<endl;
35	tim	1011
36			exit(1);
37			}
38
39			printMinizerInfo();
40
41			resetFrq = paramSet->getResetFrq();
42	tim	1015	nextResetIter = resetFrq;
43	tim	1011
44			writeFrq = paramSet->getWriteFrq();
45			nextWriteIter = writeFrq;
46
47			prevGrad = model->calcGrad();
48
49	tim	1015	direction = prevGrad;
50	tim	1011
51			maxIteration = paramSet->getMaxIteration();
52
53			for(currentIter = 0;currentIter < maxIteration; currentIter++){
54
55			// perform line search to minimize f(x + stepSize * direction) where stepSize > 0
56			lsMinimizer->minimize(direction, 0.0, 1.0);
57
58			lsStatus = lsMinimizer->getMinimizationStatus();
59
60			if(lsStatus ==MINSTATUS_ERROR){
61			minStatus = MINSTATUS_ERROR;
62			return;
63			}
64
65			prevMinX = minX;
66	tim	1015	lamda = lsMinimizer->getMinVar();
67	tim	1011
68	tim	1015	for(int i = 0; i < direction.size(); i++)
69			minX[i] = minX[i] + lamda * direction[i];
70
71	tim	1011	//calculate the gradient
72			prevGrad = gradient;
73
74			gradient = model->calcGrad();
75
76			// stop if converge
77	tim	1015	if (checkConvergence() > 0){
78	tim	1011	writeOut(minX, currentIter);
79
80			minStatus = MINSTATUS_CONVERGE;
81			return;
82			}
83
84
85			//calculate the
86	tim	1015	gamma = calcGamma(gradient, prevGrad);
87	tim	1011
88			// update new direction
89			prevDirection = direction;
90
91	tim	1015	for(int i = 0; i < direction.size(); i++)
92			direction[i] += gamma * direction[i];
93
94	tim	1011	//
95			if (currentIter == nextWriteIter){
96			nextWriteIter += writeFrq;
97			writeOut(minX, currentIter);
98			}
99
100			if (currentIter == nextResetIter){
101			reset();
102			nextResetIter += resetFrq;
103			}
104
105			}
106
107			// if writeFrq is not a multipiler of maxIteration, we need to write the final result
108			// otherwise, we already write it inside the loop, just skip it
109			if(currentIter != (nextWriteIter - writeFrq))
110			writeOut(minX, currentIter);
111
112			minStatus = MINSTATUS_MAXITER;
113			return;
114			}
115
116			int ConjugateMinimizerBase::checkConvergence(){
117
118			//test absolute gradient tolerance
119
120	tim	1015	if (sqrt(dot(gradient, gradient)) < paramSet->getGradTol())
121	tim	1011	return 1;
122			else
123			return -1;
124			}
125
126			void ConjugateMinimizerBase::reset(){
127
128			}
129
130			double FRCGMinimizer::calcGamma(vector<double>& newGrad, vector<double>& oldGrad){
131	tim	1015	return dot(newGrad, newGrad) / dot(oldGrad, newGrad);
132	tim	1011	}
133
134			double PRCGMinimizer::calcGamma(vector<double>& newGrad, vector<double>& oldGrad){
135			double gamma;
136			vector<double> deltaGrad;
137	tim	1015
138			for(int i = 0; i < newGrad.size(); i++)
139			deltaGrad.push_back(newGrad[i] - oldGrad[i]);
140	tim	1011
141	tim	1015	return dot(deltaGrad, newGrad) / dot(oldGrad, oldGrad);
142	tim	1011
143			}