OOPSE/libmdtools/NLModel1.cpp

#include "NLModel.hpp"
#include "Utility.hpp"


//calculate hessian using finite difference
SymMatrix NLModel1::FiniteHessian(vector<double>& x, vector<double>& h){

  SymMatrix hessian(ndim);
  vector<double> gplus;
  vector<double> tempX;  
  double hi;
  
#ifdef IS_MPI
  vector<double> gplusAll;
  vector<double> tempGradAll;

  //tempGradAll = getAllGrad();
#endif

  tempX = x;
    
  for(int i = 0; i < ndim; i++){    

#ifndef IS_MPI
    hi = copysign(h[i], tempX[i]);

    tempX[i] += hi;

    gplus = calcGrad(tempX);

    //hessian.Colume(i) = (gplus - currentGrad) / hi;  

    //restore tempX to its original value
    tempX[i] -= hi;
#else
    if(procMappingArray[i] == myRank){

      hi = copysign(h[i], tempX[i]);
      tempX[i] += hi;

    }
    
    gplus = calcGrad(tempX);

    gplusAll.clear();

    //very inefficient, need to change
    for(int j = 0; j , numOfProc; j++){
      //brocastVector(gplus, j);
      //gplusAll.insert(gplus);
    }
    
    if(procMappingArray[i] == myRank){

      hi = copysign(h[i], tempX[i]);
      tempX[i] += hi;

    }

    //hessian.Cloume(i) = (gplusAll - tempGradAll) / hi;  
#endif

  }

  return hessian;
  
}

//----------------------------------------------------------------------------//
ConcreteNLModel1::ConcreteNLModel1(int dim, ObjFunctor1* func ,  ConstraintList* cons)
                         : NLModel1(dim, cons){
  objfunc = func; 
}
double ConcreteNLModel1::calcF(){
  
  currentF = (*objfunc)(currentX, currentGrad);
  numOfFunEval ++;

  return currentF;
}

double ConcreteNLModel1::calcF(vector<double>& x){

  vector<double> tempGrad(x.size());
  
  double tempF;

  tempF = (*objfunc)(x, tempGrad);
  numOfFunEval ++;
  
  return tempF;  
}

vector<double> ConcreteNLModel1::calcGrad(){
  
  currentF = (*objfunc)(currentX, currentGrad);

  return currentGrad;

}

vector<double> ConcreteNLModel1::calcGrad(vector<double>& x){
  vector<double> tempGrad(x.size());

  double tempF;
  
  tempF = (*objfunc)(x, tempGrad);

  return tempGrad;
}


SymMatrix ConcreteNLModel1::calcHessian(){
  calcGrad(currentX);
  
  return FiniteHessian(currentX, currentX);
}

SymMatrix ConcreteNLModel1::calcHessian(vector<double>& x){
  return FiniteHessian(x, x);
}

Revision:	1031
Committed:	Fri Feb 6 18:58:06 2004 UTC (20 years, 5 months ago) by tim
File size:	2360 byte(s)
Log Message:	Add some lines into global.cpp to make it work with energy minimization
#	User	Rev	Content
1	tim	996	#include "NLModel.hpp"
2	tim	1023	#include "Utility.hpp"
3
4
5
6	tim	996	//calculate hessian using finite difference
7			SymMatrix NLModel1::FiniteHessian(vector<double>& x, vector<double>& h){
8
9			SymMatrix hessian(ndim);
10			vector<double> gplus;
11			vector<double> tempX;
12			double hi;
13
14			#ifdef IS_MPI
15			vector<double> gplusAll;
16	tim	1000	vector<double> tempGradAll;
17	tim	996
18	tim	1031	//tempGradAll = getAllGrad();
19	tim	996	#endif
20
21			tempX = x;
22
23			for(int i = 0; i < ndim; i++){
24
25			#ifndef IS_MPI
26			hi = copysign(h[i], tempX[i]);
27
28			tempX[i] += hi;
29
30			gplus = calcGrad(tempX);
31
32	tim	1010	//hessian.Colume(i) = (gplus - currentGrad) / hi;
33	tim	996
34			//restore tempX to its original value
35			tempX[i] -= hi;
36			#else
37			if(procMappingArray[i] == myRank){
38
39			hi = copysign(h[i], tempX[i]);
40			tempX[i] += hi;
41
42			}
43
44			gplus = calcGrad(tempX);
45
46			gplusAll.clear();
47
48			//very inefficient, need to change
49			for(int j = 0; j , numOfProc; j++){
50	tim	1031	//brocastVector(gplus, j);
51			//gplusAll.insert(gplus);
52	tim	996	}
53
54			if(procMappingArray[i] == myRank){
55
56			hi = copysign(h[i], tempX[i]);
57			tempX[i] += hi;
58
59			}
60
61	tim	1031	//hessian.Cloume(i) = (gplusAll - tempGradAll) / hi;
62	tim	996	#endif
63
64			}
65
66			return hessian;
67
68			}
69
70	tim	1023	//----------------------------------------------------------------------------//
71			ConcreteNLModel1::ConcreteNLModel1(int dim, ObjFunctor1* func , ConstraintList* cons)
72			: NLModel1(dim, cons){
73			objfunc = func;
74			}
75			double ConcreteNLModel1::calcF(){
76	tim	1000
77			currentF = (*objfunc)(currentX, currentGrad);
78	tim	1010	numOfFunEval ++;
79	tim	996
80	tim	1010	return currentF;
81	tim	996	}
82
83	tim	1023	double ConcreteNLModel1::calcF(vector<double>& x){
84	tim	996
85	tim	1000	vector<double> tempGrad(x.size());
86
87			double tempF;
88
89	tim	1015	tempF = (*objfunc)(x, tempGrad);
90	tim	1000	numOfFunEval ++;
91
92			return tempF;
93	tim	996	}
94
95	tim	1023	vector<double> ConcreteNLModel1::calcGrad(){
96	tim	1010
97	tim	1000	currentF = (*objfunc)(currentX, currentGrad);
98
99			return currentGrad;
100
101	tim	996	}
102
103	tim	1023	vector<double> ConcreteNLModel1::calcGrad(vector<double>& x){
104	tim	1015	vector<double> tempGrad(x.size());
105	tim	996
106	tim	1000	double tempF;
107
108	tim	1015	tempF = (*objfunc)(x, tempGrad);
109	tim	1000
110			return tempGrad;
111	tim	996	}
112
113	tim	1023
114			SymMatrix ConcreteNLModel1::calcHessian(){
115	tim	1000	calcGrad(currentX);
116	tim	1023
117			return FiniteHessian(currentX, currentX);
118	tim	996	}
119
120	tim	1023	SymMatrix ConcreteNLModel1::calcHessian(vector<double>& x){
121			return FiniteHessian(x, x);
122	tim	996	}
123	tim	1023