source: trunk/test_more/for_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.4 KB
Line 
1/* $Id: for_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
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
396} // End empty namespace
397
398bool for_sparse_jac(void)
399{       bool ok = true;
400
401        ok &= case_one();
402        ok &= case_two();
403        ok &= case_three();
404
405        return ok;
406}
Note: See TracBrowser for help on using the repository browser.