ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/NLModel0.cpp
Revision: 996
Committed: Wed Jan 28 22:44:44 2004 UTC (21 years, 9 months ago) by tim
File size: 4532 byte(s)
Log Message: 
Revision of NLModel0 and NLModel1

File Contents

# User Rev Content
1 tim 996 #include "NLModel.hpp"
2    
3    
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 = (*objfunc)(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] = x[i] + hi;
38     tempX[i] -= hi;
39     }
40    
41     fminus = (*objfunc)(tempX);
42    
43     if(procMappingArray[i] == myRank){
44     partialGrad[i] = (fx - fminus) / hi;
45    
46     //restore tempX to its original value
47     tempX[i] += hi;
48     }
49    
50     #endif
51     }
52    
53     return partialGrad;
54     }
55    
56     /**
57     * calculate gradient using forward finite difference
58     * df(Xk)/dXi = (f(Xk+h*ei) - f(Xk)) /h
59     * where h is a small positive scalar and ei is the ith unit vector (ith column of the identity Matrix)
60     * h can be used for all paritial derivatives, but in some cases, it is essential to use a different value
61     * for each partial derivative
62     */
63     vector<double> NLModel0::ForwardGrad(const vector<double>& x, double& fx, vector<double>& grad, const vector<double>& h){
64     vector<double> tempX = x;
65     vector<double> partialGrad;
66     double fplus;
67     double hi;
68    
69     for(int i = 0; i < ndim; i++){
70    
71     #ifndef IS_MPI
72     hi = copysign(h[i], tempX[i]);
73    
74     //tempX[i] = x[i] + hi;
75     tempX[i] += hi;
76    
77     fplus = (*objfunc)(tempX);
78    
79     partialGrad[i] = (fplus - fx) / hi;
80    
81     //restore tempX to its original value
82     tempX[i] -= hi;
83     #else
84    
85     if(procMappingArray[i] == myRank){
86    
87     hi = copysign(h[i], tempX[i]);
88    
89     //tempX[i] = x[i] + hi;
90     tempX[i] += hi;
91     }
92    
93     fminus = (*objfunc)(tempX);
94    
95     if(procMappingArray[i] == myRank){
96     partialGrad[i] = (fx - fminus) / hi;
97    
98     //restore tempX to its original value
99     tempX[i] -= hi;
100     }
101    
102     #endif
103     }
104    
105     return partialGrad;
106     }
107    
108     /**
109     * calculate gradient using central finite difference
110     * df(Xk)/dXi = (f(Xk+h*ei) - f(Xk - h*ei )) /h
111     * where h is a small positive scalar and ei is the ith unit vector (ith column of the identity Matrix)
112     * h can be used for all paritial derivatives, but in some cases, it is essential to use a different value
113     * for each partial derivative
114     */
115     vector<double> NLModel0::CentralGrad(const vector<double>& x, double& fx, vector<double>& grad, const vector<double>& h){
116     vector<double> tempX = x;
117     vector<double> partialGrad;
118     double fplus, fminus;
119     double hi;
120    
121     for(int i = 0; i < ndim; i++){
122    
123     #ifndef IS_MPI
124     hi = copysign(h[i], tempX[i]);
125    
126     //tempX[i] = x[i] + hi
127     tempX[i] += hi;
128    
129     fplus = (*objfunc)(tempX);
130    
131     //tempX[i] = x[i] -hi
132     tempX[i] -= 2*hi;
133     fminus = (*objfunc)(tempX);
134    
135     partialGrad[i] = (fplus + fminus) / (2*hi);
136    
137     //restore tempX to its original value
138     tempX[i] += hi;
139     #else
140    
141     if(procMappingArray[i] == myRank){
142    
143     hi = copysign(h[i], tempX[i]);
144    
145     //tempX[i] = x[i] + hi;
146     tempX[i] += hi;
147     }
148    
149     fplus = (*objfunc)(tempX);
150    
151     if(procMappingArray[i] == myRank){
152     partialGrad[i] = (fx - fminus) / hi;
153    
154     //restore tempX to its original value
155     tempX[i] -= 2*hi;
156     }
157    
158     fminus = (*objfunc)(tempX);
159    
160     if(procMappingArray[i] == myRank){
161     partialGrad[i] = (fx - fminus) / (2*hi);
162    
163     //restore tempX to its original value
164     tempX[i] -= hi;
165     }
166     #endif
167     }
168    
169     return partialGrad;
170     }
171    
172     /**
173     * calculate hessian using finite difference
174     * d2f(Xk)/dxidxj = (df(Xk+h*ej)/dXi + df(Xk - h*ej)/dXi) /2h
175     * where h is a small positive scalar and ei is the ith unit vector (ith column of the identity Matrix)
176     */
177     SymMatrix NLModel0::FDHessian(vector<double>& h){
178     SymMatrix H(ndim);
179    
180     for(int i = 0; i < ndim; i++){
181    
182     for(int j = i + 1; j < ndim; j++){
183    
184     }
185     }
186    
187     }

Properties

Name Value
svn:executable *