source: trunk/test_more/rev_sparse_jac.cpp @ 2506

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

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

  • Property svn:keywords set to Id
File size: 9.3 KB
Line 
1/* $Id: rev_sparse_jac.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
14# include <cppad/cppad.hpp>
15# include <vector>
16# include <valarray>
17
18# define CheckOp(Op)                   \
19        Y[index] = X[0] Op 2.;         \
20        Check[index * n + 0] = true;   \
21        Check[index * n + 1] = false;  \
22        Check[index * n + 2] = false;  \
23        index++;                       \
24        Y[index] = X[0] Op X[1];       \
25        Check[index * n + 0] = true;   \
26        Check[index * n + 1] = true;   \
27        Check[index * n + 2] = false;  \
28        index++;                       \
29        Y[index] = 3.   Op X[1];       \
30        Check[index * n + 0] = false;  \
31        Check[index * n + 1] = true;   \
32        Check[index * n + 2] = false;  \
33        index++;
34
35# define CheckUnaryFun(Fun)            \
36        Y[index] = Fun(X[0]);          \
37        Check[index * n + 0] = true;   \
38        Check[index * n + 1] = false;  \
39        Check[index * n + 2] = false;  \
40        index++;                       \
41        Y[index] = Fun(X[0] + X[1]);   \
42        Check[index * n + 0] = true;   \
43        Check[index * n + 1] = true;   \
44        Check[index * n + 2] = false;  \
45        index++;                       \
46        Y[index] = Fun(X[1]);          \
47        Check[index * n + 0] = false;  \
48        Check[index * n + 1] = true;   \
49        Check[index * n + 2] = false;  \
50        index++;
51 
52
53# define CheckBinaryFun(Fun)           \
54        Y[index] = Fun( X[0] , 2.);    \
55        Check[index * n + 0] = true;   \
56        Check[index * n + 1] = false;  \
57        Check[index * n + 2] = false;  \
58        index++;                       \
59        Y[index] = Fun( X[0] , X[1]);  \
60        Check[index * n + 0] = true;   \
61        Check[index * n + 1] = true;   \
62        Check[index * n + 2] = false;  \
63        index++;                       \
64        Y[index] = Fun( 3.   , X[1]);  \
65        Check[index * n + 0] = false;  \
66        Check[index * n + 1] = true;   \
67        Check[index * n + 2] = false;  \
68        index++;
69
70namespace { // BEGIN empty namespace
71
72bool case_one()
73{       bool ok = true;
74        using namespace CppAD;
75
76        // dimension of the domain space
77        size_t n = 3; 
78
79        // dimension of the range space
80        size_t m = (4 + 11 + 1) * 3 + 4;
81
82        // independent variable vector
83        CPPAD_TESTVECTOR(AD<double>) X(n);
84        X[0] = .1; 
85        X[1] = .2;
86        X[2] = .3;
87        Independent(X);
88
89        // dependent variable vector
90        CPPAD_TESTVECTOR(AD<double>) Y(m);
91
92        // check results vector
93        CPPAD_TESTVECTOR( bool )       Check(m * n);
94
95        // initialize index into Y
96        size_t index = 0;
97
98        // 4 binary operators
99        CheckOp(+);
100        CheckOp(-);
101        CheckOp(*);
102        CheckOp(/);
103
104        // 11 unary functions
105        CheckUnaryFun(abs);
106        CheckUnaryFun(acos);
107        CheckUnaryFun(asin);
108        CheckUnaryFun(atan);
109        CheckUnaryFun(cos);
110        CheckUnaryFun(cosh);
111        CheckUnaryFun(exp);
112        CheckUnaryFun(log);
113        CheckUnaryFun(sin);
114        CheckUnaryFun(sinh);
115        CheckUnaryFun(sqrt);
116
117        // 1 binary function
118        CheckBinaryFun(pow);
119
120        // conditional expression
121        Y[index] = CondExpLt(X[0], X[1], X[0], AD<double>(2.));
122        Check[index * n + 0] = true;
123        Check[index * n + 1] = false;
124        Check[index * n + 2] = false;
125        index++;
126        Y[index] = CondExpLt(X[0], X[1], X[0], X[1]);
127        Check[index * n + 0] = true;
128        Check[index * n + 1] = true;
129        Check[index * n + 2] = false;
130        index++;
131        Y[index] = CondExpLt(X[0], X[1], AD<double>(3.), X[1]);
132        Check[index * n + 0] = false;
133        Check[index * n + 1] = true;
134        Check[index * n + 2] = false;
135        index++;
136
137        // non-trival composition
138        Y[index] = Y[0] + Y[1] + X[2];
139        Check[index * n + 0] = true;
140        Check[index * n + 1] = true;
141        Check[index * n + 2] = true;
142        index++;
143
144        // check final index
145        assert( index == m );
146
147
148        // create function object F : X -> Y
149        ADFun<double> F(X, Y);
150
151        // --------------------------------------------------------
152        // dependency matrix for the identity function U(y) = y
153        CPPAD_TESTVECTOR( bool ) Py(m * m);
154        size_t i, j;
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
161        // evaluate the dependency matrix for F(x)
162        CPPAD_TESTVECTOR( bool ) Px(m * n);
163        Px = F.RevSparseJac(m, Py);
164
165        // check values
166        for(i = 0; i < m; i++)
167        {       for(j = 0; j < n; j++)
168                        ok &= (Px[i * n + j] == Check[i * n + j]);
169        }       
170        // --------------------------------------------------------
171        // dependency matrix for the identity function U(y) = y
172        CPPAD_TESTVECTOR(std::set<size_t>) Sy(m);
173        for(i = 0; i < m; i++)
174        {       assert( Sy[i].empty() );
175                Sy[i].insert(i);
176        }
177
178        // evaluate the dependency matrix for U(F(x))
179        CPPAD_TESTVECTOR(std::set<size_t>) Sx(m);
180        Sx = F.RevSparseJac(m, Sy);
181
182        // check values
183        std::set<size_t>::iterator itr;
184        bool found;
185        for(i = 0; i < m; i++)
186        {       for(j = 0; j < n; j++)
187                {       found = Sx[i].find(j) != Sx[i].end();
188                        ok    &= (found == Check[i * n + j]);
189                }
190        }       
191
192        return ok;
193}
194
195bool case_two()
196{       bool ok = true;
197        using namespace CppAD;
198
199        // dimension of the domain space
200        size_t n = 3; 
201
202        // dimension of the range space
203        size_t m = 3;
204
205        // inialize the vector as zero
206        CppAD::VecAD<double> Z(n - 1);
207        size_t k;
208        for(k = 0; k < n-1; k++)
209                Z[k] = 0.;
210
211        // independent variable vector
212        CPPAD_TESTVECTOR(AD<double>) X(n);
213        X[0] = 0.; 
214        X[1] = 1.;
215        X[2] = 2.;
216        Independent(X);
217
218        // VecAD vector is going to depend on X[1] and X[2]
219        Z[ X[0] ] = X[1];
220        Z[ X[1] ] = X[2]; 
221
222        // dependent variable vector
223        CPPAD_TESTVECTOR(AD<double>) Y(m);
224
225        // check results vector
226        CPPAD_TESTVECTOR( bool )       Check(m * n);
227
228        // initialize index into Y
229        size_t index = 0;
230
231        // First component only depends on X[0];
232        Y[index]             = X[0];
233        Check[index * n + 0] = true;
234        Check[index * n + 1] = false;
235        Check[index * n + 2] = false;
236        index++;
237
238        // Second component depends on the vector Z
239        AD<double> zero(0);
240        Y[index]             = Z[zero]; // Load by a parameter
241        Check[index * n + 0] = false;
242        Check[index * n + 1] = true;
243        Check[index * n + 2] = true;
244        index++;
245
246        // Third component depends on the vector Z
247        Y[index]             = Z[ X[0] ]; // Load by a variable
248        Check[index * n + 0] = false;
249        Check[index * n + 1] = true;
250        Check[index * n + 2] = true;
251        index++;
252
253        // check final index
254        assert( index == m );
255
256        // create function object F : X -> Y
257        ADFun<double> F(X, Y);
258
259        // dependency matrix for the identity function S(y) = y
260        CPPAD_TESTVECTOR( bool ) Py(m * m);
261        size_t i, j;
262        for(i = 0; i < m; i++)
263        {       for(j = 0; j < m; j++)
264                        Py[ i * m + j ] = false;
265                Py[ i * m + i ] = true;
266        }
267
268        // evaluate the dependency matrix for S [ F(x) ]
269        CPPAD_TESTVECTOR( bool ) Px(m * n);
270        Px = F.RevSparseJac(m, Py);
271
272        // check values
273        for(i = 0; i < m; i++)
274        {       for(j = 0; j < n; j++)
275                        ok &= (Px[i * n + j] == Check[i * n + j]);
276        }       
277        // --------------------------------------------------------
278        // dependency matrix for the identity function U(y) = y
279        CPPAD_TESTVECTOR(std::set<size_t>) Sy(m);
280        for(i = 0; i < m; i++)
281        {       assert( Sy[i].empty() );
282                Sy[i].insert(i);
283        }
284
285        // evaluate the dependency matrix for U(F(x))
286        CPPAD_TESTVECTOR(std::set<size_t>) Sx(m);
287        Sx = F.RevSparseJac(m, Sy);
288
289        // check values
290        std::set<size_t>::iterator itr;
291        bool found;
292        for(i = 0; i < m; i++)
293        {       for(j = 0; j < n; j++)
294                {       found = Sx[i].find(j) != Sx[i].end();
295                        ok    &= (found == Check[i * n + j]);
296                }
297        }       
298
299        return ok;
300}
301
302bool case_three()
303{       bool ok = true;
304        using namespace CppAD;
305
306        // dimension of the domain space
307        size_t n = 2; 
308
309        // dimension of the range space
310        size_t m = 3;
311
312        // independent variable vector
313        CPPAD_TESTVECTOR(AD<double>) X(n);
314        X[0] = 2.; 
315        X[1] = 3.;
316        Independent(X);
317
318        // dependent variable vector
319        CPPAD_TESTVECTOR(AD<double>) Y(m);
320
321        // check results vector
322        CPPAD_TESTVECTOR( bool )       Check(m * n);
323
324        // initialize index into Y
325        size_t index = 0;
326
327        // Y[0] only depends on X[0];
328        Y[index]             = pow(X[0], 2.);
329        Check[index * n + 0] = true;
330        Check[index * n + 1] = false;
331        index++;
332
333        // Y[1] depends on X[1]
334        Y[index]             = pow(2., X[1]);
335        Check[index * n + 0] = false;
336        Check[index * n + 1] = true;
337        index++;
338
339        // Y[2] depends on X[0] and X[1]
340        Y[index]             = pow(X[0], X[1]);
341        Check[index * n + 0] = true;
342        Check[index * n + 1] = true;
343        index++;
344
345        // check final index
346        assert( index == m );
347
348        // create function object F : X -> Y
349        ADFun<double> F(X, Y);
350
351        // -----------------------------------------------------------------
352        // dependency matrix for the identity function
353        CPPAD_TESTVECTOR( bool ) Py(m * m);
354        size_t i, j;
355        for(i = 0; i < m; i++)
356        {       for(j = 0; j < m; j++)
357                        Py[ i * m + j ] = (i == j);
358        }
359
360        // evaluate the dependency matrix for F(x)
361        CPPAD_TESTVECTOR( bool ) Px(m * n);
362        Px = F.RevSparseJac(m, Py);
363
364        // check values
365        for(i = 0; i < m; i++)
366        {       for(j = 0; j < n; j++)
367                        ok &= (Px[i * n + j] == Check[i * n + j]);
368        }       
369
370        // ---------------------------------------------------------
371        // dependency matrix for the identity function
372        CPPAD_TESTVECTOR(std::set<size_t>) Sy(m);
373        for(i = 0; i < m; i++)
374        {       assert( Sy[i].empty() );
375                Sy[i].insert(i);
376        }
377
378        // evaluate the dependency matrix for F(x)
379        CPPAD_TESTVECTOR(std::set<size_t>) Sx(m);
380        Sx = F.RevSparseJac(m, Sy);
381
382        // check values
383        bool found;
384        for(i = 0; i < m; i++)
385        {       for(j = 0; j < n; j++)
386                {       found = Sx[i].find(j) != Sx[i].end();
387                        ok &= (found == Check[i * n + j]);
388                }
389        }       
390
391        return ok;
392}
393
394} // END empty namespace
395
396bool rev_sparse_jac(void)
397{       bool ok = true;
398
399        ok &= case_one();
400        ok &= case_two();
401        ok &= case_three();
402
403        return ok;
404}
Note: See TracBrowser for help on using the repository browser.