ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/NLModel0.cpp
Revision: 1015
Committed: Tue Feb 3 22:54:52 2004 UTC (20 years, 5 months ago) by tim
File size: 6201 byte(s)
Log Message: 
NLModel0, NLModel1 pass uit test

File Contents

# User Rev Content
1 tim 996 #include "NLModel.hpp"
2 tim 1010 #include <math.h>
3 tim 996
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 tim 1010 fminus = calcF(tempX);
26 tim 996
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 tim 1010 fminus = calcF(tempX);
41 tim 996
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 tim 1010 fplus = calcF(tempX);
77 tim 996
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 tim 1015 fplus = calcF(tempX);
93 tim 996
94     if(procMappingArray[i] == myRank){
95 tim 1015 partialGrad[i] = (fplus - fx) / hi;
96 tim 996
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 tim 1010 fplus = calcF(tempX);
129 tim 996
130     //tempX[i] = x[i] -hi
131     tempX[i] -= 2*hi;
132 tim 1010 fminus = calcF(tempX);
133 tim 996
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 tim 1010 fplus = calcF(tempX);
149 tim 996
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 tim 1015 fminus = calcF(tempX);
158 tim 996
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 tim 1000 *
176 tim 996 */
177 tim 1000 SymMatrix NLModel0::FiniteHessian(vector<double>& x, double fx, vector<double>& h){
178 tim 996 SymMatrix H(ndim);
179 tim 1000 vector<double> tempX;
180     vector<double> fi(ndim);
181     vector<double> hi;
182 tim 1010 double fii, fij;
183 tim 1000
184     tempX = x;
185 tim 996
186 tim 1000 //calculate f(X+ h*ei)
187 tim 996 for(int i = 0; i < ndim; i++){
188    
189 tim 1000 #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 tim 1010 tempX[i] -= hi[i];
198 tim 1000
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 tim 1010 H(i, i) = ((fx - fi[i]) + (fii - fi[i])) / (hi[i]*hi[i]);
226 tim 1000
227     tempX[i] -= hi[i];
228    
229 tim 996 for(int j = i + 1; j < ndim; j++){
230    
231 tim 1000 tempX[j] += hi[j];
232    
233     fij = calcF(tempX);
234    
235 tim 1010 H(i,j) = ((fx - fi[i]) + (fij - fi[j])) / (hi[i]*hi[j]);
236 tim 1000
237     tempX[j] -= hi[j];
238 tim 996 }
239 tim 1000
240     tempX[i] -= hi[i];
241 tim 996 }
242 tim 1000
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 tim 1010 H(i,i) = ((fx - fi[i]) + (fii - fi[i])) / (hi[i]*hi[i]);
253 tim 1000
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 tim 1010 H(i,j) = ((fx - fi[i]) + (fij - fi[j])) / (hi[i]*hi[j]);
266 tim 1000
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 tim 996 }

Properties

Name Value
svn:executable *