source: trunk/test_more/sparse_jacobian.cpp @ 1658

Last change on this file since 1658 was 1658, checked in by bradbell, 10 years ago

/home/bradbell/cppad/trunk: Fix bool packed bug in case where set size is a multiple of n_bit.

spars_jacobian.cpp: add test for case where set size is a multiple of n_bit.
*/makefile.in: update makefile.am version corresponding to this file.
whats_new_10.omh: user's view of the changes.
sparse_pack.hpp: avoid accessing of end of data_ vector.

  • Property svn:keywords set to Id
File size: 8.9 KB
Line 
1/* $Id: sparse_jacobian.cpp 1658 2010-03-03 16:42:54Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-10 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/*
14Old sparse Jacobian example
15*/
16
17# include <cppad/cppad.hpp>
18namespace { // ---------------------------------------------------------
19
20template <class VectorBase, class VectorBool> 
21bool reverse_bool()
22{       bool ok = true;
23        using CppAD::AD;
24        using CppAD::NearEqual;
25        size_t i, j, k;
26
27        // domain space vector
28        size_t n = 4;
29        CPPAD_TEST_VECTOR< AD<double> >  X(n);
30        for(j = 0; j < n; j++)
31                X[j] = AD<double> (0);
32
33        // declare independent variables and starting recording
34        CppAD::Independent(X);
35
36        size_t m = 3;
37        CPPAD_TEST_VECTOR< AD<double> >  Y(m);
38        Y[0] = X[0] + X[1];
39        Y[1] = X[2] + X[3];
40        Y[2] = X[0] + X[1] + X[2] + X[3] * X[3] / 2.;
41
42        // create f: X -> Y and stop tape recording
43        CppAD::ADFun<double> f(X, Y);
44
45        // new value for the independent variable vector
46        VectorBase x(n);
47        for(j = 0; j < n; j++)
48                x[j] = double(j);
49
50        // Jacobian of y without sparsity pattern
51        VectorBase jac(m * n);
52        jac = f.SparseJacobian(x);
53        /*
54              [ 1 1 0 0  ]
55        jac = [ 0 0 1 1  ]
56              [ 1 1 1 x_3]
57        */
58        VectorBase check(m * n);
59        check[0] = 1.; check[1] = 1.; check[2]  = 0.; check[3]  = 0.;
60        check[4] = 0.; check[5] = 0.; check[6]  = 1.; check[7]  = 1.;
61        check[8] = 1.; check[9] = 1.; check[10] = 1.; check[11] = x[3];
62        for(k = 0; k < 12; k++)
63                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
64
65        // test passing sparsity pattern
66        VectorBool s(m * m);
67        VectorBool p(m * n);
68        for(i = 0; i < m; i++)
69        {       for(k = 0; k < m; k++)
70                        s[i * m + k] = false;
71                s[i * m + i] = true;
72        }
73        p   = f.RevSparseJac(m, s);
74        jac = f.SparseJacobian(x);
75        for(k = 0; k < 12; k++)
76                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
77
78        return ok;
79}
80
81template <class VectorBase, class VectorSet> 
82bool reverse_set()
83{       bool ok = true;
84        using CppAD::AD;
85        using CppAD::NearEqual;
86        size_t i, j, k;
87
88        // domain space vector
89        size_t n = 4;
90        CPPAD_TEST_VECTOR< AD<double> >  X(n);
91        for(j = 0; j < n; j++)
92                X[j] = AD<double> (0);
93
94        // declare independent variables and starting recording
95        CppAD::Independent(X);
96
97        size_t m = 3;
98        CPPAD_TEST_VECTOR< AD<double> >  Y(m);
99        Y[0] = X[0] + X[1];
100        Y[1] = X[2] + X[3];
101        Y[2] = X[0] + X[1] + X[2] + X[3] * X[3] / 2.;
102
103        // create f: X -> Y and stop tape recording
104        CppAD::ADFun<double> f(X, Y);
105
106        // new value for the independent variable vector
107        VectorBase x(n);
108        for(j = 0; j < n; j++)
109                x[j] = double(j);
110
111        // Jacobian of y without sparsity pattern
112        VectorBase jac(m * n);
113        jac = f.SparseJacobian(x);
114        /*
115              [ 1 1 0 0  ]
116        jac = [ 0 0 1 1  ]
117              [ 1 1 1 x_3]
118        */
119        VectorBase check(m * n);
120        check[0] = 1.; check[1] = 1.; check[2]  = 0.; check[3]  = 0.;
121        check[4] = 0.; check[5] = 0.; check[6]  = 1.; check[7]  = 1.;
122        check[8] = 1.; check[9] = 1.; check[10] = 1.; check[11] = x[3];
123        for(k = 0; k < 12; k++)
124                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
125
126        // test passing sparsity pattern
127        VectorSet s(m), p(m);
128        for(i = 0; i < m; i++)
129                s[i].insert(i);
130        p   = f.RevSparseJac(m, s);
131        jac = f.SparseJacobian(x);
132        for(k = 0; k < 12; k++)
133                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
134
135        return ok;
136}
137
138template <class VectorBase, class VectorBool> 
139bool forward_bool()
140{       bool ok = true;
141        using CppAD::AD;
142        using CppAD::NearEqual;
143        size_t j, k;
144
145        // domain space vector
146        size_t n = 3;
147        CPPAD_TEST_VECTOR< AD<double> >  X(n);
148        for(j = 0; j < n; j++)
149                X[j] = AD<double> (0);
150
151        // declare independent variables and starting recording
152        CppAD::Independent(X);
153
154        size_t m = 4;
155        CPPAD_TEST_VECTOR< AD<double> >  Y(m);
156        Y[0] = X[0] + X[2];
157        Y[1] = X[0] + X[2];
158        Y[2] = X[1] + X[2];
159        Y[3] = X[1] + X[2] * X[2] / 2.;
160
161        // create f: X -> Y and stop tape recording
162        CppAD::ADFun<double> f(X, Y);
163
164        // new value for the independent variable vector
165        VectorBase x(n);
166        for(j = 0; j < n; j++)
167                x[j] = double(j);
168
169        // Jacobian of y without sparsity pattern
170        VectorBase jac(m * n);
171        jac = f.SparseJacobian(x);
172        /*
173              [ 1 0 1   ]
174        jac = [ 1 0 1   ]
175              [ 0 1 1   ]
176              [ 0 1 x_2 ]
177        */
178        VectorBase check(m * n);
179        check[0] = 1.; check[1]  = 0.; check[2]  = 1.; 
180        check[3] = 1.; check[4]  = 0.; check[5]  = 1.;
181        check[6] = 0.; check[7]  = 1.; check[8]  = 1.; 
182        check[9] = 0.; check[10] = 1.; check[11] = x[2];
183        for(k = 0; k < 12; k++)
184                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
185
186        // test passing sparsity pattern
187        VectorBool r(n * n);
188        VectorBool p(m * n);
189        for(j = 0; j < n; j++)
190        {       for(k = 0; k < n; k++)
191                        r[j * n + k] = false;
192                r[j * n + j] = true;
193        }
194        p   = f.ForSparseJac(n, r);
195        jac = f.SparseJacobian(x);
196        for(k = 0; k < 12; k++)
197                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
198
199        return ok;
200}
201
202template <class VectorBase, class VectorSet> 
203bool forward_set()
204{       bool ok = true;
205        using CppAD::AD;
206        using CppAD::NearEqual;
207        size_t j, k;
208
209        // domain space vector
210        size_t n = 3;
211        CPPAD_TEST_VECTOR< AD<double> >  X(n);
212        for(j = 0; j < n; j++)
213                X[j] = AD<double> (0);
214
215        // declare independent variables and starting recording
216        CppAD::Independent(X);
217
218        size_t m = 4;
219        CPPAD_TEST_VECTOR< AD<double> >  Y(m);
220        Y[0] = X[0] + X[2];
221        Y[1] = X[0] + X[2];
222        Y[2] = X[1] + X[2];
223        Y[3] = X[1] + X[2] * X[2] / 2.;
224
225        // create f: X -> Y and stop tape recording
226        CppAD::ADFun<double> f(X, Y);
227
228        // new value for the independent variable vector
229        VectorBase x(n);
230        for(j = 0; j < n; j++)
231                x[j] = double(j);
232
233        // Jacobian of y without sparsity pattern
234        VectorBase jac(m * n);
235        jac = f.SparseJacobian(x);
236        /*
237              [ 1 0 1   ]
238        jac = [ 1 0 1   ]
239              [ 0 1 1   ]
240              [ 0 1 x_2 ]
241        */
242        VectorBase check(m * n);
243        check[0] = 1.; check[1]  = 0.; check[2]  = 1.; 
244        check[3] = 1.; check[4]  = 0.; check[5]  = 1.;
245        check[6] = 0.; check[7]  = 1.; check[8]  = 1.; 
246        check[9] = 0.; check[10] = 1.; check[11] = x[2];
247        for(k = 0; k < 12; k++)
248                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
249
250        // test passing sparsity pattern
251        VectorSet r(n), p(m);
252        for(j = 0; j < n; j++)
253                r[j].insert(j);
254        p   = f.ForSparseJac(n, r);
255        jac = f.SparseJacobian(x);
256        for(k = 0; k < 12; k++)
257                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
258
259        return ok;
260}
261
262bool multiple_of_n_bit()
263{       bool ok = true;
264        using CppAD::AD;
265        using CppAD::vector;
266        size_t i, j;
267
268        // should be the same as the corresponding typedef in
269        // cppad/local/sparse_pack.hpp
270        typedef size_t Pack;
271
272        // number of bits per packed value
273        size_t n_bit = std::numeric_limits<Pack>::digits;
274
275        // check case where number of variables is equal to n_bit
276        vector< AD<double> > x(n_bit);
277        vector< AD<double> > y(n_bit);
278
279        // create an AD function with domain and range dimension equal to n_bit
280        CppAD::Independent(x);
281        for(i = 0; i < n_bit; i++)
282                y[i] = x[n_bit - i - 1];
283        CppAD::ADFun<double> f(x, y);
284
285        // Jacobian sparsity patterns
286        vector<bool> r(n_bit * n_bit);
287        vector<bool> s(n_bit * n_bit);
288        for(i = 0; i < n_bit; i++)
289        {       for(j = 0; j < n_bit; j++)
290                        r[ i * n_bit + j ] = (i == j);
291        }
292        s = f.ForSparseJac(n_bit, r);
293
294        // check the result
295        for(i = 0; i < n_bit; i++)
296        {       for(j = 0; j < n_bit; j++)
297                {       if( i == n_bit - j - 1 )
298                                ok = ok & s[ i * n_bit + j ];
299                        else    ok = ok & (! s[i * n_bit + j] );
300                }
301        }
302
303        return ok;
304}
305} // End empty namespace
306# include <vector>
307# include <valarray>
308bool sparse_jacobian(void)
309{       bool ok = true;
310        // ---------------------------------------------------------------
311        // vector of bool cases
312        ok &= forward_bool< CppAD::vector<double>, CppAD::vectorBool   >();
313        ok &= reverse_bool< CppAD::vector<double>, CppAD::vector<bool> >();
314        //
315        ok &= forward_bool< std::vector<double>,   std::vector<bool>   >();
316        ok &= reverse_bool< std::vector<double>,   std::valarray<bool> >();
317        //
318        ok &= forward_bool< std::valarray<double>, CppAD::vectorBool   >();
319        ok &= reverse_bool< std::valarray<double>, CppAD::vector<bool> >();
320        // ---------------------------------------------------------------
321        // vector of set cases
322        typedef std::vector< std::set<size_t> >   std_vector_set;
323        typedef std::valarray< std::set<size_t> > std_valarray_set;
324        typedef CppAD::vector< std::set<size_t> > cppad_vector_set;
325        //
326        ok &= forward_set< CppAD::vector<double>, std_vector_set   >();
327        ok &= reverse_set< std::valarray<double>, std_vector_set   >();
328        //
329        ok &= forward_set< std::vector<double>,   cppad_vector_set >();
330        ok &= reverse_set< CppAD::vector<double>, cppad_vector_set >();
331        //
332        // According to section 26.3.2.3 of the 1998 C++ standard
333        // a const valarray does not return references to its elements.
334        // ok &= forward_set< std::valarray<double>, std_valarray_set >();
335        // ok &= reverse_set< std::valarray<double>, std_valarray_set >();
336        // ---------------------------------------------------------------
337        //
338        ok &= multiple_of_n_bit();
339        return ok;
340}
Note: See TracBrowser for help on using the repository browser.