1 #ifndef FILE_AUTODIFF
2 #define FILE_AUTODIFF
3
4 /**************************************************************************/
5 /* File: autodiff.hpp */
6 /* Author: Joachim Schoeberl */
7 /* Date: 24. Oct. 02 */
8 /**************************************************************************/
9
10 // Automatic differentiation datatype
11
12
13 /**
14 Datatype for automatic differentiation.
15 Contains function value and D derivatives. Algebraic
16 operations are overloaded by using product-rule etc. etc.
17 **/
18 template <int D, typename SCAL = double>
19 class AutoDiff
20 {
21 SCAL val;
22 SCAL dval[D];
23 public:
24
25 typedef AutoDiff<D,SCAL> TELEM;
26 typedef SCAL TSCAL;
27
28
29 /// elements are undefined
AutoDiff()30 AutoDiff () throw() { };
31 // { val = 0; for (int i = 0; i < D; i++) dval[i] = 0; } // !
32
33 /// copy constructor
AutoDiff(const AutoDiff & ad2)34 AutoDiff (const AutoDiff & ad2) throw()
35 {
36 val = ad2.val;
37 for (int i = 0; i < D; i++)
38 dval[i] = ad2.dval[i];
39 }
40
41 /// initial object with constant value
AutoDiff(SCAL aval)42 AutoDiff (SCAL aval) throw()
43 {
44 val = aval;
45 for (int i = 0; i < D; i++)
46 dval[i] = 0;
47 }
48
49 /// init object with (val, e_diffindex)
AutoDiff(SCAL aval,int diffindex)50 AutoDiff (SCAL aval, int diffindex) throw()
51 {
52 val = aval;
53 for (int i = 0; i < D; i++)
54 dval[i] = 0;
55 dval[diffindex] = 1;
56 }
57
58 /// assign constant value
operator =(SCAL aval)59 AutoDiff & operator= (SCAL aval) throw()
60 {
61 val = aval;
62 for (int i = 0; i < D; i++)
63 dval[i] = 0;
64 return *this;
65 }
66
67 /// returns value
Value() const68 SCAL Value() const throw() { return val; }
69
70 /// returns partial derivative
DValue(int i) const71 SCAL DValue (int i) const throw() { return dval[i]; }
72
73 /// access value
Value()74 SCAL & Value() throw() { return val; }
75
76 /// accesses partial derivative
DValue(int i)77 SCAL & DValue (int i) throw() { return dval[i]; }
78
79 ///
operator +=(const AutoDiff<D,SCAL> & y)80 AutoDiff<D,SCAL> & operator+= (const AutoDiff<D,SCAL> & y) throw()
81 {
82 val += y.val;
83 for (int i = 0; i < D; i++)
84 dval[i] += y.dval[i];
85 return *this;
86 }
87
88 ///
operator -=(const AutoDiff<D,SCAL> & y)89 AutoDiff<D,SCAL> & operator-= (const AutoDiff<D,SCAL> & y) throw()
90 {
91 val -= y.val;
92 for (int i = 0; i < D; i++)
93 dval[i] -= y.dval[i];
94 return *this;
95
96 }
97
98 ///
operator *=(const AutoDiff<D,SCAL> & y)99 AutoDiff<D,SCAL> & operator*= (const AutoDiff<D,SCAL> & y) throw()
100 {
101 for (int i = 0; i < D; i++)
102 {
103 // dval[i] *= y.val;
104 // dval[i] += val * y.dval[i];
105 dval[i] = dval[i] * y.val + val * y.dval[i];
106 }
107 val *= y.val;
108 return *this;
109 }
110
111 ///
operator *=(const SCAL & y)112 AutoDiff<D,SCAL> & operator*= (const SCAL & y) throw()
113 {
114 val *= y;
115 for (int i = 0; i < D; i++)
116 dval[i] *= y;
117 return *this;
118 }
119
120 ///
operator /=(const SCAL & y)121 AutoDiff<D,SCAL> & operator/= (const SCAL & y) throw()
122 {
123 SCAL iy = 1.0 / y;
124 val *= iy;
125 for (int i = 0; i < D; i++)
126 dval[i] *= iy;
127 return *this;
128 }
129
130 ///
operator ==(SCAL val2)131 bool operator== (SCAL val2) throw()
132 {
133 return val == val2;
134 }
135
136 ///
operator !=(SCAL val2)137 bool operator!= (SCAL val2) throw()
138 {
139 return val != val2;
140 }
141
142 ///
operator <(SCAL val2)143 bool operator< (SCAL val2) throw()
144 {
145 return val < val2;
146 }
147
148 ///
operator >(SCAL val2)149 bool operator> (SCAL val2) throw()
150 {
151 return val > val2;
152 }
153 };
154
155
156 //@{ AutoDiff helper functions.
157
158 /// prints AutoDiff
159 template<int D, typename SCAL>
operator <<(ostream & ost,const AutoDiff<D,SCAL> & x)160 inline ostream & operator<< (ostream & ost, const AutoDiff<D,SCAL> & x)
161 {
162 ost << x.Value() << ", D = ";
163 for (int i = 0; i < D; i++)
164 ost << x.DValue(i) << " ";
165 return ost;
166 }
167
168 /// AutoDiff plus AutoDiff
169 template<int D, typename SCAL>
operator +(const AutoDiff<D,SCAL> & x,const AutoDiff<D,SCAL> & y)170 inline AutoDiff<D,SCAL> operator+ (const AutoDiff<D,SCAL> & x, const AutoDiff<D,SCAL> & y) throw()
171 {
172 AutoDiff<D,SCAL> res;
173 res.Value () = x.Value()+y.Value();
174 // AutoDiff<D,SCAL> res(x.Value()+y.Value());
175 for (int i = 0; i < D; i++)
176 res.DValue(i) = x.DValue(i) + y.DValue(i);
177 return res;
178 }
179
180
181 /// AutoDiff minus AutoDiff
182 template<int D, typename SCAL>
operator -(const AutoDiff<D,SCAL> & x,const AutoDiff<D,SCAL> & y)183 inline AutoDiff<D,SCAL> operator- (const AutoDiff<D,SCAL> & x, const AutoDiff<D,SCAL> & y) throw()
184 {
185 AutoDiff<D,SCAL> res;
186 res.Value() = x.Value()-y.Value();
187 // AutoDiff<D,SCAL> res (x.Value()-y.Value());
188 for (int i = 0; i < D; i++)
189 res.DValue(i) = x.DValue(i) - y.DValue(i);
190 return res;
191 }
192
193 /// double plus AutoDiff
194 template<int D, typename SCAL>
operator +(double x,const AutoDiff<D,SCAL> & y)195 inline AutoDiff<D,SCAL> operator+ (double x, const AutoDiff<D,SCAL> & y) throw()
196 {
197 AutoDiff<D,SCAL> res;
198 res.Value() = x+y.Value();
199 for (int i = 0; i < D; i++)
200 res.DValue(i) = y.DValue(i);
201 return res;
202 }
203
204 /// AutoDiff plus double
205 template<int D, typename SCAL>
operator +(const AutoDiff<D,SCAL> & y,double x)206 inline AutoDiff<D,SCAL> operator+ (const AutoDiff<D,SCAL> & y, double x) throw()
207 {
208 AutoDiff<D,SCAL> res;
209 res.Value() = x+y.Value();
210 for (int i = 0; i < D; i++)
211 res.DValue(i) = y.DValue(i);
212 return res;
213 }
214
215
216 /// minus AutoDiff
217 template<int D, typename SCAL>
operator -(const AutoDiff<D,SCAL> & x)218 inline AutoDiff<D,SCAL> operator- (const AutoDiff<D,SCAL> & x) throw()
219 {
220 AutoDiff<D,SCAL> res;
221 res.Value() = -x.Value();
222 for (int i = 0; i < D; i++)
223 res.DValue(i) = -x.DValue(i);
224 return res;
225 }
226
227 /// AutoDiff minus double
228 template<int D, typename SCAL>
operator -(const AutoDiff<D,SCAL> & x,double y)229 inline AutoDiff<D,SCAL> operator- (const AutoDiff<D,SCAL> & x, double y) throw()
230 {
231 AutoDiff<D,SCAL> res;
232 res.Value() = x.Value()-y;
233 for (int i = 0; i < D; i++)
234 res.DValue(i) = x.DValue(i);
235 return res;
236 }
237
238 ///
239 template<int D, typename SCAL>
operator -(double x,const AutoDiff<D,SCAL> & y)240 inline AutoDiff<D,SCAL> operator- (double x, const AutoDiff<D,SCAL> & y) throw()
241 {
242 AutoDiff<D,SCAL> res;
243 res.Value() = x-y.Value();
244 for (int i = 0; i < D; i++)
245 res.DValue(i) = -y.DValue(i);
246 return res;
247 }
248
249
250 /// double times AutoDiff
251 template<int D, typename SCAL>
operator *(double x,const AutoDiff<D,SCAL> & y)252 inline AutoDiff<D,SCAL> operator* (double x, const AutoDiff<D,SCAL> & y) throw()
253 {
254 AutoDiff<D,SCAL> res;
255 res.Value() = x*y.Value();
256 for (int i = 0; i < D; i++)
257 res.DValue(i) = x*y.DValue(i);
258 return res;
259 }
260
261 /// AutoDiff times double
262 template<int D, typename SCAL>
operator *(const AutoDiff<D,SCAL> & y,double x)263 inline AutoDiff<D,SCAL> operator* (const AutoDiff<D,SCAL> & y, double x) throw()
264 {
265 AutoDiff<D,SCAL> res;
266 res.Value() = x*y.Value();
267 for (int i = 0; i < D; i++)
268 res.DValue(i) = x*y.DValue(i);
269 return res;
270 }
271
272 /// AutoDiff times AutoDiff
273 template<int D, typename SCAL>
operator *(const AutoDiff<D,SCAL> & x,const AutoDiff<D,SCAL> & y)274 inline AutoDiff<D,SCAL> operator* (const AutoDiff<D,SCAL> & x, const AutoDiff<D,SCAL> & y) throw()
275 {
276 AutoDiff<D,SCAL> res;
277 SCAL hx = x.Value();
278 SCAL hy = y.Value();
279
280 res.Value() = hx*hy;
281 for (int i = 0; i < D; i++)
282 res.DValue(i) = hx*y.DValue(i) + hy*x.DValue(i);
283
284 return res;
285 }
286
287 /// AutoDiff times AutoDiff
288 template<int D, typename SCAL>
sqr(const AutoDiff<D,SCAL> & x)289 inline AutoDiff<D,SCAL> sqr (const AutoDiff<D,SCAL> & x) throw()
290 {
291 AutoDiff<D,SCAL> res;
292 SCAL hx = x.Value();
293 res.Value() = hx*hx;
294 hx *= 2;
295 for (int i = 0; i < D; i++)
296 res.DValue(i) = hx*x.DValue(i);
297 return res;
298 }
299
300 /// Inverse of AutoDiff
301 template<int D, typename SCAL>
Inv(const AutoDiff<D,SCAL> & x)302 inline AutoDiff<D,SCAL> Inv (const AutoDiff<D,SCAL> & x)
303 {
304 AutoDiff<D,SCAL> res(1.0 / x.Value());
305 for (int i = 0; i < D; i++)
306 res.DValue(i) = -x.DValue(i) / (x.Value() * x.Value());
307 return res;
308 }
309
310
311 /// AutoDiff div AutoDiff
312 template<int D, typename SCAL>
operator /(const AutoDiff<D,SCAL> & x,const AutoDiff<D,SCAL> & y)313 inline AutoDiff<D,SCAL> operator/ (const AutoDiff<D,SCAL> & x, const AutoDiff<D,SCAL> & y)
314 {
315 return x * Inv (y);
316 }
317
318 /// AutoDiff div double
319 template<int D, typename SCAL>
operator /(const AutoDiff<D,SCAL> & x,double y)320 inline AutoDiff<D,SCAL> operator/ (const AutoDiff<D,SCAL> & x, double y)
321 {
322 return (1/y) * x;
323 }
324
325 /// double div AutoDiff
326 template<int D, typename SCAL>
operator /(double x,const AutoDiff<D,SCAL> & y)327 inline AutoDiff<D,SCAL> operator/ (double x, const AutoDiff<D,SCAL> & y)
328 {
329 return x * Inv(y);
330 }
331
332
333
334
335 template<int D, typename SCAL>
fabs(const AutoDiff<D,SCAL> & x)336 inline AutoDiff<D,SCAL> fabs (const AutoDiff<D,SCAL> & x)
337 {
338 double abs = fabs (x.Value());
339 AutoDiff<D,SCAL> res( abs );
340 if (abs != 0.0)
341 for (int i = 0; i < D; i++)
342 res.DValue(i) = x.DValue(i) / abs;
343 else
344 for (int i = 0; i < D; i++)
345 res.DValue(i) = 0.0;
346 return res;
347 }
348
349 //@}
350
351 #endif
352