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

    model->setX(minX);
    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] = -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 != (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:	1046
Committed:	Tue Feb 10 21:33:44 2004 UTC (21 years, 2 months ago) by tim
File size:	3776 byte(s)
Log Message:	single version of energy minimization is 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
75			if(lsStatus ==MINSTATUS_ERROR){
76			minStatus = MINSTATUS_ERROR;
77			return;
78			}
79
80			prevMinX = minX;
81	tim	1015	lamda = lsMinimizer->getMinVar();
82	tim	1011
83	tim	1015	for(int i = 0; i < direction.size(); i++)
84			minX[i] = minX[i] + lamda * direction[i];
85
86	tim	1011	//calculate the gradient
87			prevGrad = gradient;
88	tim	1023
89			model->setX(minX);
90	tim	1011	gradient = model->calcGrad();
91
92			// stop if converge
93	tim	1015	if (checkConvergence() > 0){
94	tim	1011	writeOut(minX, currentIter);
95
96			minStatus = MINSTATUS_CONVERGE;
97			return;
98			}
99
100
101			//calculate the
102	tim	1015	gamma = calcGamma(gradient, prevGrad);
103	tim	1011
104			// update new direction
105			prevDirection = direction;
106
107	tim	1015	for(int i = 0; i < direction.size(); i++)
108	tim	1023	direction[i] = -gradient[i] + gamma * direction[i];
109	tim	1015
110	tim	1011	//
111			if (currentIter == nextWriteIter){
112			nextWriteIter += writeFrq;
113			writeOut(minX, currentIter);
114			}
115
116			if (currentIter == nextResetIter){
117			reset();
118			nextResetIter += resetFrq;
119			}
120
121			}
122
123			// if writeFrq is not a multipiler of maxIteration, we need to write the final result
124			// otherwise, we already write it inside the loop, just skip it
125			if(currentIter != (nextWriteIter - writeFrq))
126			writeOut(minX, currentIter);
127
128			minStatus = MINSTATUS_MAXITER;
129			return;
130			}
131
132			int ConjugateMinimizerBase::checkConvergence(){
133
134			//test absolute gradient tolerance
135
136	tim	1031	if (sqrt(dotProduct(gradient, gradient)) < paramSet->getGTol())
137	tim	1011	return 1;
138			else
139			return -1;
140			}
141
142			void ConjugateMinimizerBase::reset(){
143
144			}
145
146			double FRCGMinimizer::calcGamma(vector<double>& newGrad, vector<double>& oldGrad){
147	tim	1023	return dotProduct(newGrad, newGrad) / dotProduct(oldGrad, newGrad);
148	tim	1011	}
149
150			double PRCGMinimizer::calcGamma(vector<double>& newGrad, vector<double>& oldGrad){
151			double gamma;
152			vector<double> deltaGrad;
153	tim	1015
154			for(int i = 0; i < newGrad.size(); i++)
155			deltaGrad.push_back(newGrad[i] - oldGrad[i]);
156	tim	1011
157	tim	1023	return dotProduct(deltaGrad, newGrad) / dotProduct(oldGrad, oldGrad);
158	tim	1011
159			}