ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/NLModel0.cpp
Revision: 1023
Committed: Wed Feb 4 22:26:00 2004 UTC (20 years, 5 months ago) by tim
File size: 6223 byte(s)
Log Message: 
Fix a bunch of bugs   :-)
Single version of conjugate gradient with golden search linesearch pass a couple of
functions test. Brent's  algorithm is still broken

File Contents

# Content
1 #include "NLModel.hpp"
2 #include <math.h>
3 #include "Utility.hpp"
4 /**
5 * calculate gradient using backward finite difference
6 * df(Xk)/dXi = (f(Xk) - f(Xk - h*ei)) /h
7 * where h is a small positive scalar and ei is the ith unit vector (ith column of the identity Matrix)
8 * h can be used for all paritial derivatives, but in some cases, it is essential to use a different value
9 * for each partial derivative
10 */
11 vector<double> NLModel0::BackwardGrad(const vector<double>& x, double& fx, vector<double>& grad, const vector<double>& h){
12 vector<double> tempX = x;
13 vector<double> partialGrad;
14 double fminus;
15 double hi;
16
17 for(int i = 0; i < ndim; i++){
18
19 #ifndef IS_MPI
20 hi = copysign(h[i], tempX[i]);
21
22 //tempX[i] = x[i] + hi;
23 tempX[i] -= hi;
24
25 fminus = calcF(tempX);
26
27 partialGrad[i] = (fx - fminus) / hi;
28
29 //restore tempX to its original value
30 tempX[i] += hi;
31 #else
32
33 if(procMappingArray[i] == myRank){
34
35 hi = copysign(h[i], tempX[i]);
36
37 tempX[i] -= hi;
38 }
39
40 fminus = calcF(tempX);
41
42 if(procMappingArray[i] == myRank){
43 partialGrad[i] = (fx - fminus) / hi;
44
45 //restore tempX to its original value
46 tempX[i] += hi;
47 }
48
49 #endif
50 }
51
52 return partialGrad;
53 }
54
55 /**
56 * calculate gradient using forward finite difference
57 * df(Xk)/dXi = (f(Xk+h*ei) - f(Xk)) /h
58 * where h is a small positive scalar and ei is the ith unit vector (ith column of the identity Matrix)
59 * h can be used for all paritial derivatives, but in some cases, it is essential to use a different value
60 * for each partial derivative
61 */
62 vector<double> NLModel0::ForwardGrad(const vector<double>& x, double& fx, vector<double>& grad, const vector<double>& h){
63 vector<double> tempX = x;
64 vector<double> partialGrad;
65 double fplus;
66 double hi;
67
68 for(int i = 0; i < ndim; i++){
69
70 #ifndef IS_MPI
71 hi = copysign(h[i], tempX[i]);
72
73 //tempX[i] = x[i] + hi;
74 tempX[i] += hi;
75
76 fplus = calcF(tempX);
77
78 partialGrad[i] = (fplus - fx) / hi;
79
80 //restore tempX to its original value
81 tempX[i] -= hi;
82 #else
83
84 if(procMappingArray[i] == myRank){
85
86 hi = copysign(h[i], tempX[i]);
87
88 //tempX[i] = x[i] + hi;
89 tempX[i] += hi;
90 }
91
92 fplus = calcF(tempX);
93
94 if(procMappingArray[i] == myRank){
95 partialGrad[i] = (fplus - fx) / hi;
96
97 //restore tempX to its original value
98 tempX[i] -= hi;
99 }
100
101 #endif
102 }
103
104 return partialGrad;
105 }
106
107 /**
108 * calculate gradient using central finite difference
109 * df(Xk)/dXi = (f(Xk+h*ei) - f(Xk - h*ei )) /h
110 * where h is a small positive scalar and ei is the ith unit vector (ith column of the identity Matrix)
111 * h can be used for all paritial derivatives, but in some cases, it is essential to use a different value
112 * for each partial derivative
113 */
114 vector<double> NLModel0::CentralGrad(const vector<double>& x, double& fx, vector<double>& grad, const vector<double>& h){
115 vector<double> tempX = x;
116 vector<double> partialGrad;
117 double fplus, fminus;
118 double hi;
119
120 for(int i = 0; i < ndim; i++){
121
122 #ifndef IS_MPI
123 hi = copysign(h[i], tempX[i]);
124
125 //tempX[i] = x[i] + hi
126 tempX[i] += hi;
127
128 fplus = calcF(tempX);
129
130 //tempX[i] = x[i] -hi
131 tempX[i] -= 2*hi;
132 fminus = calcF(tempX);
133
134 partialGrad[i] = (fplus + fminus) / (2*hi);
135
136 //restore tempX to its original value
137 tempX[i] += hi;
138 #else
139
140 if(procMappingArray[i] == myRank){
141
142 hi = copysign(h[i], tempX[i]);
143
144 //tempX[i] = x[i] + hi;
145 tempX[i] += hi;
146 }
147
148 fplus = calcF(tempX);
149
150 if(procMappingArray[i] == myRank){
151 partialGrad[i] = (fx - fminus) / hi;
152
153 //restore tempX to its original value
154 tempX[i] -= 2*hi;
155 }
156
157 fminus = calcF(tempX);
158
159 if(procMappingArray[i] == myRank){
160 partialGrad[i] = (fx - fminus) / (2*hi);
161
162 //restore tempX to its original value
163 tempX[i] -= hi;
164 }
165 #endif
166 }
167
168 return partialGrad;
169 }
170
171 /**
172 * calculate hessian using finite difference
173 * d2f(Xk)/dxidxj = (df(Xk+h*ej)/dXi + df(Xk - h*ej)/dXi) /2h
174 * where h is a small positive scalar and ei is the ith unit vector (ith column of the identity Matrix)
175 *
176 */
177 SymMatrix NLModel0::FiniteHessian(vector<double>& x, double fx, vector<double>& h){
178 SymMatrix H(ndim);
179 vector<double> tempX;
180 vector<double> fi(ndim);
181 vector<double> hi;
182 double fii, fij;
183
184 tempX = x;
185
186 //calculate f(X+ h*ei)
187 for(int i = 0; i < ndim; i++){
188
189 #ifndef IS_MPI
190
191 hi[i] = copysign(h[i],tempX[i]);
192
193 tempX[i] += hi[i];
194
195 fi[i] = calcF(tempX);
196
197 tempX[i] -= hi[i];
198
199 #else
200
201 if(procMappingArray[i] == myRank){
202
203 hi[i] = copysign(h[i],tempX[i]);
204
205 tempX[i] += hi[i];
206 }
207
208 fi[i] = calcF(tempX);
209
210 if(procMappingArray[i] == myRank){
211 tempX[i] -= hi[i];
212 }
213
214 #endif
215 }
216
217 #ifndef IS_MPI
218
219 for(int i = 0; i < ndim; i++){
220
221 tempX[i] += hi[i] * 2;
222
223 fii = calcF(tempX);
224
225 H(i, i) = ((fx - fi[i]) + (fii - fi[i])) / (hi[i]*hi[i]);
226
227 tempX[i] -= hi[i];
228
229 for(int j = i + 1; j < ndim; j++){
230
231 tempX[j] += hi[j];
232
233 fij = calcF(tempX);
234
235 H(i,j) = ((fx - fi[i]) + (fij - fi[j])) / (hi[i]*hi[j]);
236
237 tempX[j] -= hi[j];
238 }
239
240 tempX[i] -= hi[i];
241 }
242
243 #else
244 for(int i = 0; i < ndim; i++){
245
246 if(procMappingArray[i] == myRank){
247 tempX[i] += hi[i] * 2;
248 }
249
250 fii = calcF(tempX);
251
252 H(i,i) = ((fx - fi[i]) + (fii - fi[i])) / (hi[i]*hi[i]);
253
254 if(procMappingArray[i] == myRank){
255 tempX[i] -= hi[i];
256 }
257
258 for(int j = i + 1; j < ndim; j++){
259
260 if(procMappingArray[j] == myRank){
261 tempX[j] += hi[j];
262 }
263
264 fij = calcF(tempX);
265 H(i,j) = ((fx - fi[i]) + (fij - fi[j])) / (hi[i]*hi[j]);
266
267 if(procMappingArray[j] == myRank){
268 tempX[j] -= hi[j];
269 }
270
271 }
272
273 if(procMappingArray[i] == myRank){
274 tempX[i] -= hi[i];
275 }
276
277 }
278
279 #endif
280
281 return H;
282 }

Properties

Name Value
svn:executable *