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 - 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:	1057
Committed:	Tue Feb 17 19:23:44 2004 UTC (21 years, 9 months ago) by tim
File size:	3780 byte(s)
Log Message:	adding function shakeF in order to remove the constraint force along bond direction
#	Content
1	#include "ConjugateMinimizer.hpp"
2	#include "Utility.hpp"
3	ConjugateMinimizerBase::ConjugateMinimizerBase(NLModel1* nlmodel, MinimizerParameterSet* param)
4	: MinimizerUsingLineSearch(param){
5	int dim;
6
7	model = nlmodel;
8	//set the dimension
9
10	#ifndef IS_MPI
11	dim = model->getDim();
12	#else
13	dim = model->getDim();
14	#endif
15	prevGrad.resize(dim);
16	gradient.resize(dim);
17	prevDirection.resize(dim);
18	direction.resize(dim);
19	}
20
21	bool ConjugateMinimizerBase::isSolvable(){
22
23	//conjuage gradient can only solve unconstrained nonlinear model
24
25	if (!model->hasConstraints())
26	return true;
27	else
28	return false;
29	}
30
31	void ConjugateMinimizerBase::printMinizerInfo(){
32
33	}
34
35	void ConjugateMinimizerBase::minimize(){
36	int maxIteration;
37	int nextResetIter;
38	int resetFrq;
39	int nextWriteIter;
40	int writeFrq;
41	int lsStatus;
42	double gamma;
43	double lamda;
44
45
46	if (!isSolvable()){
47	cout << "ConjugateMinimizerBase Error: This nonlinear model can not be solved by " << minimizerName <<endl;
48
49	exit(1);
50	}
51
52	printMinizerInfo();
53
54	resetFrq = paramSet->getResetFrq();
55	nextResetIter = resetFrq;
56
57	writeFrq = paramSet->getWriteFrq();
58	nextWriteIter = writeFrq;
59
60	minX = model->getX();
61	gradient = model->calcGrad();
62
63	for(int i = 0; i < direction.size(); i++)
64	direction[i] = -gradient[i];
65
66	maxIteration = paramSet->getMaxIteration();
67
68	for(currentIter = 1;currentIter <= maxIteration; currentIter++){
69
70	// perform line search to minimize f(x + lamda * direction) where stepSize > 0
71	lsMinimizer->minimize(direction, 0.0, 0.01);
72
73	lsStatus = lsMinimizer->getMinimizationStatus();
74
75	if(lsStatus ==MINSTATUS_ERROR){
76	minStatus = MINSTATUS_ERROR;
77	return;
78	}
79
80	prevMinX = minX;
81	lamda = lsMinimizer->getMinVar();
82
83	for(int i = 0; i < direction.size(); i++)
84	minX[i] = minX[i] + lamda * direction[i];
85
86	//calculate the gradient
87	prevGrad = gradient;
88
89	model->setX(minX);
90	gradient = model->calcGrad();
91
92	// stop if converge
93	if (checkConvergence() > 0){
94	writeOut(minX, currentIter);
95
96	minStatus = MINSTATUS_CONVERGE;
97	return;
98	}
99
100
101	//calculate the
102	gamma = calcGamma(gradient, prevGrad);
103
104	// update new direction
105	prevDirection = direction;
106
107	for(int i = 0; i < direction.size(); i++)
108	direction[i] = -gradient[i] + gamma * direction[i];
109
110	//
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 - 1 != (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	if (sqrt(dotProduct(gradient, gradient)) < paramSet->getGTol())
137	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	return dotProduct(newGrad, newGrad) / dotProduct(oldGrad, newGrad);
148	}
149
150	double PRCGMinimizer::calcGamma(vector<double>& newGrad, vector<double>& oldGrad){
151	double gamma;
152	vector<double> deltaGrad;
153
154	for(int i = 0; i < newGrad.size(); i++)
155	deltaGrad.push_back(newGrad[i] - oldGrad[i]);
156
157	return dotProduct(deltaGrad, newGrad) / dotProduct(oldGrad, oldGrad);
158
159	}