source: trunk/test_more/cond_exp_ad.cpp @ 3520

Last change on this file since 3520 was 2506, checked in by bradbell, 7 years ago

Change Licenses: CPL-1.0 -> EPL-1.0, GPL-2.0->GPL-3.0

  • Property svn:keywords set to Id
File size: 7.3 KB
Line 
1/* $Id: cond_exp_ad.cpp 2506 2012-10-24 19:36:49Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
4
5CppAD is distributed under multiple licenses. This distribution is under
6the terms of the
7                    Eclipse 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/*
14Test of CondExp with AD< AD< Base > > types
15*/
16
17# include <cppad/cppad.hpp>
18
19typedef CppAD::AD< double >     ADdouble;
20typedef CppAD::AD< ADdouble > ADADdouble;
21
22namespace { // BEGIN empty namespace
23
24bool CondExpADOne(void)
25{       bool ok = true;
26
27        using namespace CppAD;
28        size_t n = 3;
29        size_t m = 8;
30
31        // ADdouble independent variable vector
32        CPPAD_TESTVECTOR( ADdouble ) Xa(n);
33        Xa[0] = -1.;
34        Xa[1] =  0.;
35        Xa[2] =  1.;
36        Independent(Xa);
37
38        // ADdouble independent variable vector
39        CPPAD_TESTVECTOR( ADADdouble ) Xaa(n);
40        Xaa[0] = Xa[0];
41        Xaa[1] = Xa[1];
42        Xaa[2] = Xa[2];
43        Independent(Xaa);
44
45        // ADADdouble parameter
46        ADADdouble p = ADADdouble(Xa[0]);
47        ADADdouble q = ADADdouble(Xa[1]);
48        ADADdouble r = ADADdouble(Xa[2]);
49
50        // ADADdouble dependent variable vector
51        CPPAD_TESTVECTOR( ADADdouble ) Yaa(m);
52
53        // CondExp(parameter, parameter, parameter)
54        Yaa[0] = CondExp(p, q, r);
55
56        // CondExp(parameter, parameter, variable)
57        Yaa[1] = CondExp(p, q, Xaa[2]);
58
59        // CondExp(parameter, varaible, parameter)
60        Yaa[2] = CondExp(p, Xaa[1], r);
61
62        // CondExp(parameter, variable, variable)
63        Yaa[3] = CondExp(p, Xaa[1], Xaa[2]);
64
65        // CondExp(variable, variable, variable)
66        Yaa[5] = CondExp(Xaa[0], Xaa[1], Xaa[2]);
67
68        // CondExp(variable, variable, parameter)
69        Yaa[4] = CondExp(Xaa[0], Xaa[1], r);
70
71        // CondExp(variable, parameter, variable)
72        Yaa[6] =  CondExp(Xaa[0], q, Xaa[2]);
73
74        // CondExp(variable, parameter, parameter)
75        Yaa[7] =  CondExp(Xaa[0], q, r);
76
77        // create fa: Xaa -> Yaa function object
78        ADFun< ADdouble > fa(Xaa, Yaa);
79
80        // function values
81        CPPAD_TESTVECTOR( ADdouble ) Ya(m);
82        Ya  = fa.Forward(0, Xa);
83
84        // create f: Xa -> Ya function object
85        ADFun<double> f(Xa, Ya);
86
87        // check result of function evaluation
88        CPPAD_TESTVECTOR(double) x(n);
89        CPPAD_TESTVECTOR(double) y(m);
90        x[0] = 1.;
91        x[1] = 0.;
92        x[2] = -1.;
93        y = f.Forward(0, x);
94        size_t i;
95        for(i = 0; i < m; i++)
96        {       // y[i] = CondExp(x[0], x[1], x[2])
97                if( x[0] > 0 )
98                        ok &= (y[i] == x[1]);
99                else    ok &= (y[i] == x[2]);
100        }
101
102        // check forward mode derivatives
103        CPPAD_TESTVECTOR(double) dx(n);
104        CPPAD_TESTVECTOR(double) dy(m);
105        dx[0] = 1.;
106        dx[1] = 2.;
107        dx[2] = 3.;
108        dy    = f.Forward(1, dx);
109        for(i = 0; i < m; i++)
110        {       if( x[0] > 0. )
111                        ok &= (dy[i] == dx[1]);
112                else    ok &= (dy[i] == dx[2]);
113        }
114
115        // calculate Jacobian
116        CPPAD_TESTVECTOR(double) J(m * n);
117        size_t j;
118        for(i = 0; i < m; i++)
119        {       for(j = 0; j < n; j++)
120                        J[i * n + j] = 0.;
121                if( x[0] > 0. )
122                        J[i * n + 1] = 1.;
123                else    J[i * n + 2] = 1.;
124        }
125
126        // check reverse mode derivatives
127        for(i = 0; i < m; i++)
128                dy[i] = double(i);
129        dx    = f.Reverse(1, dy);
130        double sum;
131        for(j = 0; j < n; j++)
132        {       sum = 0;
133                for(i = 0; i < m; i++)
134                        sum += dy[i] * J[i * n + j];
135                ok &= (sum == dx[j]);
136        }
137
138        // forward mode computation of sparsity pattern
139        CPPAD_TESTVECTOR(bool) Px(n * n);
140        for(i = 0; i < n; i++)
141        {       for(j = 0; j < n; j++)
142                        Px[i * n + j] = false;
143                Px[i * n + i] = true;
144        }
145        CPPAD_TESTVECTOR(bool) Py(m * n);
146        Py = f.ForSparseJac(n, Px);
147        for(i = 0; i < m; i++)
148        {       ok &= Py[ i * n + 0 ] == false;
149                ok &= Py[ i * n + 1 ] == true;
150                ok &= Py[ i * n + 2 ] == true;
151        }
152
153        // reverse mode computation of sparsity pattern
154        Py.resize(m * m);
155        for(i = 0; i < m; i++)
156        {       for(j = 0; j < m; j++)
157                        Py[i * m + j] = false;
158                Py[i * m + i] = true;
159        }
160        Px.resize(m * n);
161        Px = f.RevSparseJac(m, Py);
162        for(i = 0; i < m; i++)
163        {       for(j = 0; j < n; j++)
164                        ok &= ( Px[i * n + j] == ( j > 0 ) );
165        }
166
167        return ok;
168}
169bool CondExpADTwo(void)
170{       bool ok = true;
171
172        using namespace CppAD;
173        size_t n = 3;
174        size_t m = 8;
175
176        // ADdouble independent variable vector
177        CPPAD_TESTVECTOR( ADdouble ) Xa(n);
178        Xa[0] = -1.;
179        Xa[1] =  0.;
180        Xa[2] =  1.;
181        Independent(Xa);
182
183        // use VecAD so that sparsity results are local
184        VecAD<double> Va(1);
185        ADdouble zero = 0.;
186        Va[zero]      = Xa[0];
187
188        // ADdouble independent variable vector
189        CPPAD_TESTVECTOR( ADADdouble ) Xaa(n);
190        Xaa[0] = ADdouble( Va[zero] );
191        Xaa[1] = Xa[1];
192        Xaa[2] = Xa[2];
193        Independent(Xaa);
194
195        // ADADdouble parameter
196        ADADdouble p = ADADdouble(Xa[0]);
197        ADADdouble q = ADADdouble(Xa[1]);
198        ADADdouble r = ADADdouble(Xa[2]);
199
200        // ADADdouble dependent variable vector
201        CPPAD_TESTVECTOR( ADADdouble ) Yaa(m);
202
203        // CondExp(parameter, parameter, parameter)
204        Yaa[0] = CondExp(p, q, r);
205
206        // CondExp(parameter, parameter, variable)
207        Yaa[1] = CondExp(p, q, Xaa[2]);
208
209        // CondExp(parameter, varaible, parameter)
210        Yaa[2] = CondExp(p, Xaa[1], r);
211
212        // CondExp(parameter, variable, variable)
213        Yaa[3] = CondExp(p, Xaa[1], Xaa[2]);
214
215        // CondExp(variable, variable, variable)
216        Yaa[5] = CondExp(Xaa[0], Xaa[1], Xaa[2]);
217
218        // CondExp(variable, variable, parameter)
219        Yaa[4] = CondExp(Xaa[0], Xaa[1], r);
220
221        // CondExp(variable, parameter, variable)
222        Yaa[6] =  CondExp(Xaa[0], q, Xaa[2]);
223
224        // CondExp(variable, parameter, parameter)
225        Yaa[7] =  CondExp(Xaa[0], q, r);
226
227        // create fa: Xaa -> Yaa function object
228        ADFun< ADdouble > fa(Xaa, Yaa);
229
230        // function values
231        CPPAD_TESTVECTOR( ADdouble ) Ya(m);
232        Ya  = fa.Forward(0, Xa);
233
234        // create f: Xa -> Ya function object
235        ADFun<double> f(Xa, Ya);
236
237        // check use_VecAD
238        ok &= f.use_VecAD();
239
240        // check result of function evaluation
241        CPPAD_TESTVECTOR(double) x(n);
242        CPPAD_TESTVECTOR(double) y(m);
243        x[0] = 1.;
244        x[1] = 0.;
245        x[2] = -1.;
246        y = f.Forward(0, x);
247        size_t i;
248        for(i = 0; i < m; i++)
249        {       // y[i] = CondExp(x[0], x[1], x[2])
250                if( x[0] > 0 )
251                        ok &= (y[i] == x[1]);
252                else    ok &= (y[i] == x[2]);
253        }
254
255        // check forward mode derivatives
256        CPPAD_TESTVECTOR(double) dx(n);
257        CPPAD_TESTVECTOR(double) dy(m);
258        dx[0] = 1.;
259        dx[1] = 2.;
260        dx[2] = 3.;
261        dy    = f.Forward(1, dx);
262        for(i = 0; i < m; i++)
263        {       if( x[0] > 0. )
264                        ok &= (dy[i] == dx[1]);
265                else    ok &= (dy[i] == dx[2]);
266        }
267
268        // calculate Jacobian
269        CPPAD_TESTVECTOR(double) J(m * n);
270        size_t j;
271        for(i = 0; i < m; i++)
272        {       for(j = 0; j < n; j++)
273                        J[i * n + j] = 0.;
274                if( x[0] > 0. )
275                        J[i * n + 1] = 1.;
276                else    J[i * n + 2] = 1.;
277        }
278
279        // check reverse mode derivatives
280        for(i = 0; i < m; i++)
281                dy[i] = double(i);
282        dx    = f.Reverse(1, dy);
283        double sum;
284        for(j = 0; j < n; j++)
285        {       sum = 0;
286                for(i = 0; i < m; i++)
287                        sum += dy[i] * J[i * n + j];
288                ok &= (sum == dx[j]);
289        }
290
291        // forward mode computation of sparsity pattern
292        CPPAD_TESTVECTOR(bool) Px(n * n);
293        for(i = 0; i < n; i++)
294        {       for(j = 0; j < n; j++)
295                        Px[i * n + j] = false;
296                Px[i * n + i] = true;
297        }
298        CPPAD_TESTVECTOR(bool) Py(m * n);
299        Py = f.ForSparseJac(n, Px);
300        for(i = 0; i < m; i++)
301        {       for(j = 0; j < n; j++)
302                        // sparsity pattern works for both true and false cases.
303                        ok &= ( Py[i * n + j] == (j > 0) );
304        }
305
306        // reverse mode computation of sparsity pattern
307        Py.resize(m * m);
308        for(i = 0; i < m; i++)
309        {       for(j = 0; j < m; j++)
310                        Py[i * m + j] = false;
311                Py[i * m + i] = true;
312        }
313        Px.resize(m * n);
314        Px = f.RevSparseJac(m, Py);
315        for(i = 0; i < m; i++)
316        {       for(j = 0; j < n; j++)
317                        ok &= ( Px[i * n + j] == (j > 0) );
318        }
319
320        return ok;
321}
322
323} // END empty namespace
324
325bool CondExpAD(void)
326{       bool ok = true;
327        ok     &= CondExpADOne();
328        ok     &= CondExpADTwo();
329        return ok;
330}
Note: See TracBrowser for help on using the repository browser.