OOPSE/libmdtools/ConjugateMinimizer.cpp

#include "ConjugateMinimizer.hpp"
#include "Utility.hpp"
ConjugateMinimizerBase::ConjugateMinimizerBase(NLModel1* nlmodel, MinimizerParameterSet* param)
                                : MinimizerUsingLineSearch(param){
  int dim;
  
  model = nlmodel;
  //set the dimension
  
#ifndef IS_MPI
  dim = model->getDim();
#else
  dim = model->getDim();
#endif
  prevGrad.resize(dim);    
  gradient.resize(dim);
  prevDirection.resize(dim);
  direction.resize(dim);
}

bool ConjugateMinimizerBase::isSolvable(){

  //conjuage gradient can only solve unconstrained nonlinear model

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

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;

  minX = model->getX();
  gradient = model->calcGrad();

  for(int i = 0; i < direction.size(); i++)
    direction[i] = -gradient[i];
  
  maxIteration = paramSet->getMaxIteration();

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

    // perform line search to minimize f(x + lamda * direction) where stepSize > 0
    lsMinimizer->minimize(direction, 0.0, 0.01);
    
    lsStatus = lsMinimizer->getMinimizationStatus();

    lamda = lsMinimizer->getMinVar();
    
    if(lsStatus ==MINSTATUS_ERROR){
      if (lamda == 0){

          for(int i = 0; i < direction.size(); i++)
          direction[i] = -prevGrad[i]; 
          
        continue;
      }
        minStatus = MINSTATUS_ERROR;
      return;
    }
    else{
      prevMinX = minX;
    }

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

    //calculate the gradient
    prevGrad = gradient;

    model->setX(minX);
    gradient = model->calcGrad();

    minX = model->getX();
    // 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] = -gradient[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 - 1 != (nextWriteIter - writeFrq))
    writeOut(minX, currentIter);

  minStatus = MINSTATUS_MAXITER;
  return;
}

int ConjugateMinimizerBase::checkConvergence(){

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

void ConjugateMinimizerBase::reset(){  
  
}

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