source: trunk/test_more/cond_exp.cpp @ 1601

Last change on this file since 1601 was 1447, checked in by bradbell, 11 years ago

trunk: Merge in branches/sweep from revision 1404 to 1446.

  • Property svn:keywords set to Id
File size: 11.8 KB
Line 
1/* $Id: cond_exp.cpp 1447 2009-07-04 19:15:09Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 Bradley M. Bell
4
5CppAD is distributed under multiple licenses. This distribution is under
6the terms of the
7                    Common Public License Version 1.0.
8
9A copy of this license is included in the COPYING file of this distribution.
10Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
11-------------------------------------------------------------------------- */
12
13/*
14Comprehensive test built on 08/07 for new user interface to CondExp
15*/
16// BEGIN PROGRAM
17
18# include <cppad/cppad.hpp>
19
20namespace { // Begin empty namespace
21
22bool CondExp_pvvv(void)
23{       bool ok = true;
24
25        using namespace CppAD;
26
27        // independent variable vector
28        CPPAD_TEST_VECTOR< AD<double> > X(3);
29        X[0]     = 0.;
30        X[1]     = 1.;
31        X[2]     = 2.;
32        Independent(X);
33
34        // parameter value
35        AD<double> one = 1.; 
36
37        // dependent variable vector
38        CPPAD_TEST_VECTOR< AD<double> > Y(5);
39
40        // CondExp(parameter, variable, variable, variable)
41        Y[0] = CondExpLt(one, X[0], X[1], X[2]);
42        Y[1] = CondExpLe(one, X[0], X[1], X[2]);
43        Y[2] = CondExpEq(one, X[0], X[1], X[2]);
44        Y[3] = CondExpGe(one, X[0], X[1], X[2]);
45        Y[4] = CondExpGt(one, X[0], X[1], X[2]);
46
47        // create f: X -> Y
48        ADFun<double> f(X, Y);
49
50        // vectors for function values
51        CPPAD_TEST_VECTOR<double> v( f.Domain() );
52        CPPAD_TEST_VECTOR<double> w( f.Range() );
53
54        // vectors for derivative values
55        CPPAD_TEST_VECTOR<double> dv( f.Domain() );
56        CPPAD_TEST_VECTOR<double> dw( f.Range() );
57
58        // check original function values
59        ok &= Y[0] == X[2];
60        ok &= Y[1] == X[2];
61        ok &= Y[2] == X[2];
62        ok &= Y[3] == X[1];
63        ok &= Y[4] == X[1];
64
65        // function values
66        v[0] = 2.;
67        v[1] = 1.;
68        v[2] = 0.;
69        w    = f.Forward(0, v);
70        ok &= ( w[0] == v[1] );
71        ok &= ( w[1] == v[1] );
72        ok &= ( w[2] == v[2] );
73        ok &= ( w[3] == v[2] );
74        ok &= ( w[4] == v[2] );
75
76        // forward mode derivative values
77        dv[0] = 1.;
78        dv[1] = 2.;
79        dv[2] = 3.;
80        dw    = f.Forward(1, dv);
81        ok   &= (dw[0] == dv[1] );
82        ok   &= (dw[1] == dv[1] );
83        ok   &= (dw[2] == dv[2] );
84        ok   &= (dw[3] == dv[2] );
85        ok   &= (dw[4] == dv[2] );
86
87        // reverse mode derivative values
88        dw[0] = 1.;
89        dw[1] = 2.;
90        dw[2] = 3.;
91        dw[3] = 4.;
92        dw[4] = 5.;
93        dv    = f.Reverse(1, dw);
94        ok   &= (dv[0] == 0.);
95        ok   &= (dv[1] == dw[0] + dw[1] );
96        ok   &= (dv[2] == dw[2] + dw[3] + dw[4] );
97       
98        return ok;
99}
100bool CondExp_vpvv(void)
101{       bool ok = true;
102
103        using namespace CppAD;
104
105        // independent variable vector
106        CPPAD_TEST_VECTOR< AD<double> > X(3);
107        X[0]     = 0.;
108        X[1]     = 1.;
109        X[2]     = 2.;
110        Independent(X);
111
112        // parameter value
113        AD<double> one = 1.; 
114
115        // dependent variable vector
116        CPPAD_TEST_VECTOR< AD<double> > Y(5);
117
118        // CondExp(variable, parameter, variable, variable)
119        Y[0] = CondExpLt(X[0], one, X[1], X[2]);
120        Y[1] = CondExpLe(X[0], one, X[1], X[2]);
121        Y[2] = CondExpEq(X[0], one, X[1], X[2]);
122        Y[3] = CondExpGe(X[0], one, X[1], X[2]);
123        Y[4] = CondExpGt(X[0], one, X[1], X[2]);
124
125        // create f: X -> Y
126        ADFun<double> f(X, Y);
127
128        // vectors for function values
129        CPPAD_TEST_VECTOR<double> v( f.Domain() );
130        CPPAD_TEST_VECTOR<double> w( f.Range() );
131
132        // vectors for derivative values
133        CPPAD_TEST_VECTOR<double> dv( f.Domain() );
134        CPPAD_TEST_VECTOR<double> dw( f.Range() );
135
136        // check original function values
137        ok &= Y[0] == X[1];
138        ok &= Y[1] == X[1];
139        ok &= Y[2] == X[2];
140        ok &= Y[3] == X[2];
141        ok &= Y[4] == X[2];
142
143        // function values
144        v[0] = 2.;
145        v[1] = 1.;
146        v[2] = 0.;
147        w    = f.Forward(0, v);
148        ok &= ( w[0] == v[2] );
149        ok &= ( w[1] == v[2] );
150        ok &= ( w[2] == v[2] );
151        ok &= ( w[3] == v[1] );
152        ok &= ( w[4] == v[1] );
153
154        // forward mode derivative values
155        dv[0] = 1.;
156        dv[1] = 2.;
157        dv[2] = 3.;
158        dw    = f.Forward(1, dv);
159        ok   &= (dw[0] == dv[2] );
160        ok   &= (dw[1] == dv[2] );
161        ok   &= (dw[2] == dv[2] );
162        ok   &= (dw[3] == dv[1] );
163        ok   &= (dw[4] == dv[1] );
164
165        // reverse mode derivative values
166        dw[0] = 1.;
167        dw[1] = 2.;
168        dw[2] = 3.;
169        dw[3] = 4.;
170        dw[4] = 5.;
171        dv    = f.Reverse(1, dw);
172        ok   &= (dv[0] == 0.);
173        ok   &= (dv[1] == dw[3] + dw[4] );
174        ok   &= (dv[2] == dw[0] + dw[1] + dw[2] );
175       
176        return ok;
177}
178bool CondExp_vvpv(void)
179{       bool ok = true;
180
181        using namespace CppAD;
182
183        // independent variable vector
184        CPPAD_TEST_VECTOR< AD<double> > X(3);
185        X[0]     = 0.;
186        X[1]     = 1.;
187        X[2]     = 2.;
188        Independent(X);
189
190        // parameter value
191        AD<double> three = 3.; 
192
193        // dependent variable vector
194        CPPAD_TEST_VECTOR< AD<double> > Y(5);
195
196        // CondExp(variable, variable, parameter, variable)
197        Y[0] = CondExpLt(X[0], X[1], three, X[2]);
198        Y[1] = CondExpLe(X[0], X[1], three, X[2]);
199        Y[2] = CondExpEq(X[0], X[1], three, X[2]);
200        Y[3] = CondExpGe(X[0], X[1], three, X[2]);
201        Y[4] = CondExpGt(X[0], X[1], three, X[2]);
202
203        // create f: X -> Y
204        ADFun<double> f(X, Y);
205
206        // vectors for function values
207        CPPAD_TEST_VECTOR<double> v( f.Domain() );
208        CPPAD_TEST_VECTOR<double> w( f.Range() );
209
210        // vectors for derivative values
211        CPPAD_TEST_VECTOR<double> dv( f.Domain() );
212        CPPAD_TEST_VECTOR<double> dw( f.Range() );
213
214        // check original function values
215        ok &= Y[0] == three;
216        ok &= Y[1] == three;
217        ok &= Y[2] == X[2];
218        ok &= Y[3] == X[2];
219        ok &= Y[4] == X[2];
220
221        // function values
222        v[0] = 2.;
223        v[1] = 1.;
224        v[2] = 0.;
225        w    = f.Forward(0, v);
226        ok &= ( w[0] == v[2] );
227        ok &= ( w[1] == v[2] );
228        ok &= ( w[2] == v[2] );
229        ok &= ( w[3] == three );
230        ok &= ( w[4] == three );
231
232        // forward mode derivative values
233        dv[0] = 1.;
234        dv[1] = 2.;
235        dv[2] = 3.;
236        dw    = f.Forward(1, dv);
237        ok   &= (dw[0] == dv[2] );
238        ok   &= (dw[1] == dv[2] );
239        ok   &= (dw[2] == dv[2] );
240        ok   &= (dw[3] == 0.    );
241        ok   &= (dw[4] == 0.    );
242
243        // reverse mode derivative values
244        dw[0] = 1.;
245        dw[1] = 2.;
246        dw[2] = 3.;
247        dw[3] = 4.;
248        dw[4] = 5.;
249        dv    = f.Reverse(1, dw);
250        ok   &= (dv[0] == 0.);
251        ok   &= (dv[1] == 0.);
252        ok   &= (dv[2] == dw[0] + dw[1] + dw[2] );
253       
254        return ok;
255}
256bool CondExp_vvvp(void)
257{       bool ok = true;
258
259        using namespace CppAD;
260
261        // independent variable vector
262        CPPAD_TEST_VECTOR< AD<double> > X(3);
263        X[0]     = 0.;
264        X[1]     = 1.;
265        X[2]     = 2.;
266        Independent(X);
267
268        // parameter value
269        AD<double> three = 3.; 
270
271        // dependent variable vector
272        CPPAD_TEST_VECTOR< AD<double> > Y(5);
273
274        // CondExp(variable, variable, variable, parameter)
275        Y[0] = CondExpLt(X[0], X[1], X[2], three);
276        Y[1] = CondExpLe(X[0], X[1], X[2], three);
277        Y[2] = CondExpEq(X[0], X[1], X[2], three);
278        Y[3] = CondExpGe(X[0], X[1], X[2], three);
279        Y[4] = CondExpGt(X[0], X[1], X[2], three);
280
281        // create f: X -> Y
282        ADFun<double> f(X, Y);
283
284        // vectors for function values
285        CPPAD_TEST_VECTOR<double> v( f.Domain() );
286        CPPAD_TEST_VECTOR<double> w( f.Range() );
287
288        // vectors for derivative values
289        CPPAD_TEST_VECTOR<double> dv( f.Domain() );
290        CPPAD_TEST_VECTOR<double> dw( f.Range() );
291
292        // check original function values
293        ok &= Y[0] == X[2];
294        ok &= Y[1] == X[2];
295        ok &= Y[2] == three;
296        ok &= Y[3] == three;
297        ok &= Y[4] == three;
298
299        // function values
300        v[0] = 2.;
301        v[1] = 1.;
302        v[2] = 0.;
303        w    = f.Forward(0, v);
304        ok &= ( w[0] == three );
305        ok &= ( w[1] == three );
306        ok &= ( w[2] == three );
307        ok &= ( w[3] == v[2]  );
308        ok &= ( w[4] == v[2]  );
309
310        // forward mode derivative values
311        dv[0] = 1.;
312        dv[1] = 2.;
313        dv[2] = 3.;
314        dw    = f.Forward(1, dv);
315        ok   &= (dw[0] == 0.    );
316        ok   &= (dw[1] == 0.    );
317        ok   &= (dw[2] == 0.    );
318        ok   &= (dw[3] == dv[2] );
319        ok   &= (dw[4] == dv[2] );
320
321        // reverse mode derivative values
322        dw[0] = 1.;
323        dw[1] = 2.;
324        dw[2] = 3.;
325        dw[3] = 4.;
326        dw[4] = 5.;
327        dv    = f.Reverse(1, dw);
328        ok   &= (dv[0] == 0.);
329        ok   &= (dv[1] == 0.);
330        ok   &= (dv[2] == dw[3] + dw[4] );
331       
332        return ok;
333}
334
335# include <limits>
336bool SecondOrderReverse(void)
337{       // Bradley M. Bell 2009-07-04
338        // Reverse mode for CExpOp was only modifying the highest order partial
339        // This test demonstrated the bug
340        bool ok = true;
341        using CppAD::AD;
342        using CppAD::NearEqual;
343        double eps = 10. * std::numeric_limits<double>::epsilon();
344
345        size_t n = 1;
346        CPPAD_TEST_VECTOR< AD<double> > X(n);
347        X[0] = 2.;
348        CppAD::Independent(X);
349
350        size_t m = 2;
351        CPPAD_TEST_VECTOR< AD<double> > Y(m);
352
353        AD<double> left = X[0];
354        AD<double> right = X[0] * X[0]; 
355        AD<double> exp_if_true  = left;
356        AD<double> exp_if_false = right;
357
358        // Test of reverse mode using exp_if_true case
359        // For this value of X, should be the same as Z = X[0]
360        AD<double> Z = CondExpLt(left, right, exp_if_true, exp_if_false);
361        Y[0] = Z * Z;
362
363        // Test of reverse mode  using exp_if_false case
364        exp_if_false = left;
365        exp_if_true  = right;
366        Z            = CondExpGt(left, right, exp_if_true, exp_if_false);
367        Y[1]         = Z * Z;
368       
369        CppAD::ADFun<double> f(X, Y);
370
371        // first order forward
372        CPPAD_TEST_VECTOR<double> dx(n);
373        size_t p = 1;
374        dx[0]    = 1.;
375        f.Forward(p, dx);
376
377        // second order reverse (test exp_if_true case)
378        CPPAD_TEST_VECTOR<double> w(m), dw(2 * n);
379        w[0] = 1.;
380        w[1] = 0.;
381        p    = 2;
382        dw = f.Reverse(p, w);
383
384        // check first derivative in dw
385        double check = 2. * Value( X[0] );
386        ok &= NearEqual(dw[0], check, eps, eps); 
387
388        // check second derivative in dw
389        check = 2.;
390        ok &= NearEqual(dw[1], check, eps, eps); 
391
392        // test exp_if_false case
393        w[0] = 0.;
394        w[1] = 1.;
395        p    = 2;
396        dw = f.Reverse(p, w);
397
398        // check first derivative in dw
399        check = 2. * Value( X[0] );
400        ok &= NearEqual(dw[0], check, eps, eps); 
401
402        // check second derivative in dw
403        check = 2.;
404        ok &= NearEqual(dw[1], check, eps, eps); 
405
406        return ok;
407}
408
409double Infinity(double zero)
410{       return 1. / zero; }
411
412bool OldExample(void)
413{       bool ok = true;
414
415        using CppAD::AD;
416        using CppAD::NearEqual;
417        using CppAD::log; 
418        using CppAD::abs;
419        double eps = 100. * std::numeric_limits<double>::epsilon();
420
421        // domain space vector
422        size_t n = 5;
423        CPPAD_TEST_VECTOR< AD<double> > X(n);
424        size_t j;
425        for(j = 0; j < n; j++)
426                X[j] = 1.;
427
428        // declare independent variables and start tape recording
429        CppAD::Independent(X);
430
431        // sum with respect to j of log of absolute value of X[j]
432        // sould be - infinity if any of the X[j] are zero
433        AD<double> MinusInfinity = - Infinity(0.);
434        AD<double> Sum           = 0.;
435        AD<double> Zero(0);
436        for(j = 0; j < n; j++)
437        {       // if X[j] > 0
438                Sum += CppAD::CondExpGt(X[j], Zero, log(X[j]),     Zero);
439
440                // if X[j] < 0
441                Sum += CppAD::CondExpLt(X[j], Zero, log(-X[j]),    Zero);
442
443                // if X[j] == 0
444                Sum += CppAD::CondExpEq(X[j], Zero, MinusInfinity, Zero);
445        }
446
447        // range space vector
448        size_t m = 1;
449        CPPAD_TEST_VECTOR< AD<double> > Y(m);
450        Y[0] = Sum;
451
452        // create f: X -> Y and stop tape recording
453        CppAD::ADFun<double> f(X, Y);
454
455        // vectors for arguments to the function object f
456        CPPAD_TEST_VECTOR<double> x(n);   // argument values
457        CPPAD_TEST_VECTOR<double> y(m);   // function values
458        CPPAD_TEST_VECTOR<double> w(m);   // function weights
459        CPPAD_TEST_VECTOR<double> dw(n);  // derivative of weighted function
460
461        // a case where abs( x[j] ) > 0 for all j
462        double check  = 0.;
463        double sign   = 1.;
464        for(j = 0; j < n; j++)
465        {       sign *= -1.;
466                x[j] = sign * double(j + 1); 
467                check += log( abs( x[j] ) );
468        }
469
470        // function value
471        y  = f.Forward(0, x);
472        ok &= ( y[0] == check );
473
474        // compute derivative of y[0]
475        w[0] = 1.;
476        dw   = f.Reverse(1, w);
477        for(j = 0; j < n; j++)
478        {       if( x[j] > 0. )
479                        ok &= NearEqual(dw[j], 1./abs( x[j] ), eps, eps); 
480                else    ok &= NearEqual(dw[j], -1./abs( x[j] ), eps, eps); 
481        }
482
483        // a case where x[0] is equal to zero
484        sign = 1.;
485        for(j = 0; j < n; j++)
486        {       sign *= -1.;
487                x[j] = sign * double(j); 
488        }
489
490        // function value
491        y   = f.Forward(0, x);
492        ok &= ( y[0] == -Infinity(0.) );
493
494        // compute derivative of y[0]
495        w[0] = 1.;
496        dw   = f.Reverse(1, w);
497        for(j = 0; j < n; j++)
498        {       if( x[j] > 0. )
499                        ok &= NearEqual(dw[j], 1./abs( x[j] ), eps, eps); 
500                else if( x[j] < 0. )
501                        ok &= NearEqual(dw[j], -1./abs( x[j] ), eps, eps); 
502                else
503                {       // in this case computing dw[j] ends up multiplying
504                        // -infinity * zero and hence results in Nan
505                }
506        }
507       
508        return ok;
509}
510} // end empty namespace
511
512bool CondExp(void)
513{       bool ok  = true;
514        ok      &= CondExp_pvvv();
515        ok      &= CondExp_vpvv();
516        ok      &= CondExp_vvpv();
517        ok      &= CondExp_vvvp();
518        ok      &= SecondOrderReverse();
519        ok      &= OldExample();
520        return ok;
521}
522// END PROGRAM
Note: See TracBrowser for help on using the repository browser.