source: trunk/test_more/for_sparse_jac.cpp @ 3520

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

merge in changes from branches/atomic; see bin/svn_merge.sh

  • Property svn:keywords set to Id
File size: 11.4 KB
Line 
1/* $Id: for_sparse_jac.cpp 2859 2013-05-28 06:03:21Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 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
70
71namespace { // Begin empty namespace
72
73bool case_one()
74{       bool ok = true;
75        using namespace CppAD;
76
77        // dimension of the domain space
78        size_t n = 3; 
79
80        // dimension of the range space
81        size_t m = (4 + 11 + 1) * 3 + 4;
82
83        // independent variable vector
84        CPPAD_TESTVECTOR(AD<double>) X(n);
85        X[0] = .1; 
86        X[1] = .2;
87        X[2] = .3;
88        Independent(X);
89
90        // dependent variable vector
91        CPPAD_TESTVECTOR(AD<double>) Y(m);
92
93        // check results vector
94        CPPAD_TESTVECTOR( bool )       Check(m * n);
95
96        // initialize index into Y
97        size_t index = 0;
98
99        // 4 binary operators
100        CheckOp(+);
101        CheckOp(-);
102        CheckOp(*);
103        CheckOp(/);
104
105        // 11 unary functions
106        CheckUnaryFun(abs);
107        CheckUnaryFun(acos);
108        CheckUnaryFun(asin);
109        CheckUnaryFun(atan);
110        CheckUnaryFun(cos);
111        CheckUnaryFun(cosh);
112        CheckUnaryFun(exp);
113        CheckUnaryFun(log);
114        CheckUnaryFun(sin);
115        CheckUnaryFun(sinh);
116        CheckUnaryFun(sqrt);
117
118        // 1 binary function
119        CheckBinaryFun(pow);
120
121        // conditional expression (value of comparision does not matter)
122        Y[index] = CondExpLt(X[0], X[1], X[0], AD<double>(2.));
123        Check[index * n + 0] = true;
124        Check[index * n + 1] = false;
125        Check[index * n + 2] = false;
126        index++;
127        Y[index] = CondExpLt(X[0], X[1], X[0], X[1]);
128        Check[index * n + 0] = true;
129        Check[index * n + 1] = true;
130        Check[index * n + 2] = false;
131        index++;
132        Y[index] = CondExpLt(X[0], X[1], AD<double>(3.), X[1]);
133        Check[index * n + 0] = false;
134        Check[index * n + 1] = true;
135        Check[index * n + 2] = false;
136        index++;
137
138        // non-trival composition
139        Y[index] = X[0] * X[1] + X[1] * X[2];
140        Check[index * n + 0] = true;
141        Check[index * n + 1] = true;
142        Check[index * n + 2] = true;
143        index++;
144
145        // check final index
146        assert( index == m );
147
148        // create function object F : X -> Y
149        ADFun<double> F(X, Y);
150
151        // ---------------------------------------------------------
152        // dependency matrix for the identity function W(x) = x
153        CPPAD_TESTVECTOR( bool ) Px(n * n);
154        size_t i, j;
155        for(i = 0; i < n; i++)
156        {       for(j = 0; j < n; j++)
157                        Px[ i * n + j ] = false;
158                Px[ i * n + i ] = true;
159        }
160
161        // evaluate the dependency matrix for F(X(x))
162        CPPAD_TESTVECTOR( bool ) Py(m * n);
163        Py = F.ForSparseJac(n, Px);
164
165        // check values
166        for(i = 0; i < m; i++)
167        {       for(j = 0; j < n; j++)
168                        ok &= (Py[i * n + j] == Check[i * n + j]);
169        }       
170
171        // ---------------------------------------------------------
172        // dependency matrix for the identity function W(x) = x
173        CPPAD_TESTVECTOR(std::set<size_t>) Sx(n);
174        for(i = 0; i < n; i++)
175        {       assert( Sx[i].empty() );
176                Sx[i].insert(i);
177        }
178
179        // evaluate the dependency matrix for F(X(x))
180        CPPAD_TESTVECTOR(std::set<size_t>) Sy(m);
181        Sy = F.ForSparseJac(n, Sx);
182
183        // check values
184        bool found;
185        for(i = 0; i < m; i++)
186        {       for(j = 0; j < n; j++)
187                {       found = Sy[i].find(j) != Sy[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        // -----------------------------------------------------------------
260        // dependency matrix for the identity function W(x) = x
261        CPPAD_TESTVECTOR( bool ) Px(n * n);
262        size_t i, j;
263        for(i = 0; i < n; i++)
264        {       for(j = 0; j < n; j++)
265                        Px[ i * n + j ] = false;
266                Px[ i * n + i ] = true;
267        }
268
269        // evaluate the dependency matrix for F(X(x))
270        CPPAD_TESTVECTOR( bool ) Py(m * n);
271        Py = F.ForSparseJac(n, Px);
272
273        // check values
274        for(i = 0; i < m; i++)
275        {       for(j = 0; j < n; j++)
276                        ok &= (Py[i * n + j] == Check[i * n + j]);
277        }       
278
279        // ---------------------------------------------------------
280        // dependency matrix for the identity function W(x) = x
281        CPPAD_TESTVECTOR(std::set<size_t>) Sx(n);
282        for(i = 0; i < n; i++)
283        {       assert( Sx[i].empty() );
284                Sx[i].insert(i);
285        }
286
287        // evaluate the dependency matrix for F(X(x))
288        CPPAD_TESTVECTOR(std::set<size_t>) Sy(m);
289        Sy = F.ForSparseJac(n, Sx);
290
291        // check values
292        bool found;
293        for(i = 0; i < m; i++)
294        {       for(j = 0; j < n; j++)
295                {       found = Sy[i].find(j) != Sy[i].end();
296                        ok &= (found == Check[i * n + j]);
297                }
298        }       
299
300        return ok;
301}
302
303bool case_three()
304{       bool ok = true;
305        using namespace CppAD;
306
307        // dimension of the domain space
308        size_t n = 2; 
309
310        // dimension of the range space
311        size_t m = 3;
312
313        // independent variable vector
314        CPPAD_TESTVECTOR(AD<double>) X(n);
315        X[0] = 2.; 
316        X[1] = 3.;
317        Independent(X);
318
319        // dependent variable vector
320        CPPAD_TESTVECTOR(AD<double>) Y(m);
321
322        // check results vector
323        CPPAD_TESTVECTOR( bool )       Check(m * n);
324
325        // initialize index into Y
326        size_t index = 0;
327
328        // Y[0] only depends on X[0];
329        Y[index]             = pow(X[0], 2.);
330        Check[index * n + 0] = true;
331        Check[index * n + 1] = false;
332        index++;
333
334        // Y[1] depends on X[1]
335        Y[index]             = pow(2., X[1]);
336        Check[index * n + 0] = false;
337        Check[index * n + 1] = true;
338        index++;
339
340        // Y[2] depends on X[0] and X[1]
341        Y[index]             = pow(X[0], X[1]);
342        Check[index * n + 0] = true;
343        Check[index * n + 1] = true;
344        index++;
345
346        // check final index
347        assert( index == m );
348
349        // create function object F : X -> Y
350        ADFun<double> F(X, Y);
351
352        // -----------------------------------------------------------------
353        // dependency matrix for the identity function
354        CPPAD_TESTVECTOR( bool ) Px(n * n);
355        size_t i, j;
356        for(i = 0; i < n; i++)
357        {       for(j = 0; j < n; j++)
358                        Px[ i * n + j ] = false;
359                Px[ i * n + i ] = true;
360        }
361
362        // evaluate the dependency matrix for F(X(x))
363        CPPAD_TESTVECTOR( bool ) Py(m * n);
364        Py = F.ForSparseJac(n, Px);
365
366        // check values
367        for(i = 0; i < m; i++)
368        {       for(j = 0; j < n; j++)
369                        ok &= (Py[i * n + j] == Check[i * n + j]);
370        }       
371
372        // ---------------------------------------------------------
373        // dependency matrix for the identity function
374        CPPAD_TESTVECTOR(std::set<size_t>) Sx(n);
375        for(i = 0; i < n; i++)
376        {       assert( Sx[i].empty() );
377                Sx[i].insert(i);
378        }
379
380        // evaluate the dependency matrix for F(X(x))
381        CPPAD_TESTVECTOR(std::set<size_t>) Sy(m);
382        Sy = F.ForSparseJac(n, Sx);
383
384        // check values
385        bool found;
386        for(i = 0; i < m; i++)
387        {       for(j = 0; j < n; j++)
388                {       found = Sy[i].find(j) != Sy[i].end();
389                        ok &= (found == Check[i * n + j]);
390                }
391        }       
392
393        return ok;
394}
395
396bool case_four()
397{       bool ok = true;
398        using namespace CppAD;
399
400        // dimension of the domain space
401        size_t n = 2; 
402
403        // dimension of the range space
404        size_t m = 3;
405
406        // independent variable vector
407        CPPAD_TESTVECTOR(AD<double>) X(n);
408        X[0] = 2.; 
409        X[1] = 3.;
410        Independent(X);
411
412        // dependent variable vector
413        CPPAD_TESTVECTOR(AD<double>) Y(m);
414
415        // check results vector
416        CPPAD_TESTVECTOR( bool )       Check(m * n);
417
418        // initialize index into Y
419        size_t index = 0;
420
421        // Y[0] only depends on X[0];
422        Y[index]             = pow(X[0], 2.);
423        Check[index * n + 0] = true;
424        Check[index * n + 1] = false;
425        index++;
426
427        // Y[1] depends on X[1]
428        Y[index]             = pow(2., X[1]);
429        Check[index * n + 0] = false;
430        Check[index * n + 1] = true;
431        index++;
432
433        // Y[2] depends on X[0] and X[1]
434        Y[index]             = pow(X[0], X[1]);
435        Check[index * n + 0] = true;
436        Check[index * n + 1] = true;
437        index++;
438
439        // check final index
440        assert( index == m );
441
442        // create function object F : X -> Y
443        ADFun<double> F(X, Y);
444
445        // -----------------------------------------------------------------
446        // dependency matrix for the identity function
447        CPPAD_TESTVECTOR( bool ) Px(n * n);
448        size_t i, j;
449        for(i = 0; i < n; i++)
450        {       for(j = 0; j < n; j++)
451                        Px[ i * n + j ] = false;
452                Px[ i * n + i ] = true;
453        }
454
455        // evaluate the dependency matrix for F(X(x))
456        bool transpose = true;
457        CPPAD_TESTVECTOR( bool ) Py(n * m);
458        Py = F.ForSparseJac(n, Px, transpose);
459
460        // check values
461        for(i = 0; i < m; i++)
462        {       for(j = 0; j < n; j++)
463                        ok &= (Py[j * m + i] == Check[i * n + j]);
464        }       
465
466        // ---------------------------------------------------------
467        // dependency matrix for the identity function
468        CPPAD_TESTVECTOR(std::set<size_t>) Sx(n);
469        for(i = 0; i < n; i++)
470        {       assert( Sx[i].empty() );
471                Sx[i].insert(i);
472        }
473
474        // evaluate the dependency matrix for F(X(x))
475        CPPAD_TESTVECTOR(std::set<size_t>) Sy(n);
476        Sy = F.ForSparseJac(n, Sx, transpose);
477
478        // check values
479        bool found;
480        for(i = 0; i < m; i++)
481        {       for(j = 0; j < n; j++)
482                {       found = Sy[j].find(i) != Sy[j].end();
483                        ok &= (found == Check[i * n + j]);
484                }
485        }       
486
487        return ok;
488}
489
490} // End empty namespace
491
492bool for_sparse_jac(void)
493{       bool ok = true;
494
495        ok &= case_one();
496        ok &= case_two();
497        ok &= case_three();
498        ok &= case_four();
499
500        return ok;
501}
Note: See TracBrowser for help on using the repository browser.