source: trunk/test_more/sparse_jacobian.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: 15.0 KB
Line 
1/* $Id: sparse_jacobian.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/*
14Old sparse Jacobian example
15*/
16
17# include <cppad/cppad.hpp>
18namespace { // ---------------------------------------------------------
19
20bool rc_tridiagonal(void)
21{       bool ok = true;
22        using CppAD::AD;
23        using CppAD::NearEqual;
24        size_t i, j, k, ell;
25        double eps = 10. * CppAD::epsilon<double>();
26
27        // domain space vector
28        size_t n = 13; // must be greater than or equal 3 (see n_sweep below)
29        CPPAD_TESTVECTOR(AD<double>)  X(n);
30        CPPAD_TESTVECTOR(double)        x(n);
31        for(j = 0; j < n; j++)
32                X[j] = x[j] = double(j+1);
33
34        // declare independent variables and starting recording
35        CppAD::Independent(X);
36
37        size_t m = n;
38        CPPAD_TESTVECTOR(AD<double>)  Y(m);
39        CPPAD_TESTVECTOR(double) check(m * n );
40        for(ell = 0; ell < m * n; ell++)
41                check[ell] = 0.0;
42
43        size_t K = 0;
44        for(i = 0; i < n; i++)
45        {       ell        = i * n + i;
46                Y[i]       = (ell+1) * 0.5 * X[i] * X[i];
47                check[ell] = (ell+1) * x[i];
48                K++;
49                if( i < n-1 )
50                {       j          = i + 1;
51                        ell        = i * n + j;
52                        Y[i]      += (ell+1) * 0.5 * X[i+1] * X[i+1];
53                        check[ell] = (ell+1) * x[i+1];
54                        K++;
55                }
56                if(i > 0 )
57                {       j          = i - 1;
58                        ell        = i * n + j;
59                        Y[i]      += (ell+1) * 0.5 * X[i-1] * X[i-1];
60                        check[ell] = (ell+1) * x[i-1];
61                }
62        }
63
64        // create f: X -> Y and stop tape recording
65        CppAD::ADFun<double> f(X, Y);
66
67        // sparsity pattern
68        CppAD::vector< std::set<size_t> > s(m), p(m);
69        for(i = 0; i < m; i++)
70                s[i].insert(i);
71        p   = f.RevSparseJac(m, s);
72
73        // Request the upper triangle of the array
74        CPPAD_TESTVECTOR(size_t) r(K), c(K);
75        CPPAD_TESTVECTOR(double) jac(K);
76        k = 0;
77        for(i = 0; i < n; i++)
78        {       r[k] = i;
79                c[k] = i;
80                k++;   
81                if( i < n-1 )
82                {       r[k] = i;
83                        c[k] = i+1;
84                        k++;
85                } 
86        }
87        ok &= K == k;
88
89        CppAD::sparse_jacobian_work work;
90        size_t n_sweep = f.SparseJacobianForward(x, p, r, c, jac, work);
91        ok &= n_sweep == 3;
92        for(k = 0; k < K; k++)
93        {       ell = r[k] * n + c[k];
94                ok &=  NearEqual(check[ell], jac[k], eps, eps);
95        }
96        work.clear();
97        n_sweep = f.SparseJacobianReverse(x, p, r, c, jac, work);
98        ok &= n_sweep == 3;
99        for(k = 0; k < K; k++)
100        {       ell = r[k] * n + c[k];
101                ok &=  NearEqual(check[ell], jac[k], eps, eps);
102        }
103
104        return ok;
105}
106
107template <class VectorBase, class VectorSet> 
108bool rc_set(void)
109{       bool ok = true;
110        using CppAD::AD;
111        using CppAD::NearEqual;
112        size_t i, j, k, ell;
113        double eps = 10. * CppAD::epsilon<double>();
114
115        // domain space vector
116        size_t n = 4;
117        CPPAD_TESTVECTOR(AD<double>)  X(n);
118        for(j = 0; j < n; j++)
119                X[j] = AD<double> (0);
120
121        // declare independent variables and starting recording
122        CppAD::Independent(X);
123
124        size_t m = 3;
125        CPPAD_TESTVECTOR(AD<double>)  Y(m);
126        Y[0] = 1.0*X[0] + 2.0*X[1];
127        Y[1] = 3.0*X[2] + 4.0*X[3];
128        Y[2] = 5.0*X[0] + 6.0*X[1] + 7.0*X[3]*X[3]/2.;
129
130        // create f: X -> Y and stop tape recording
131        CppAD::ADFun<double> f(X, Y);
132
133        // new value for the independent variable vector
134        VectorBase x(n);
135        for(j = 0; j < n; j++)
136                x[j] = double(j);
137
138        // Jacobian of y
139        /*
140              [ 1 2 0 0    ]
141        jac = [ 0 0 3 4    ]
142              [ 5 6 0 7*x_3]
143        */
144        VectorBase check(m * n);
145        check[0] = 1.; check[1] = 2.; check[2]  = 0.; check[3]  = 0.;
146        check[4] = 0.; check[5] = 0.; check[6]  = 3.; check[7]  = 4.;
147        check[8] = 5.; check[9] = 6.; check[10] = 0.; check[11] = 7.*x[3];
148
149        // sparsity pattern
150        VectorSet s(m), p(m);
151        for(i = 0; i < m; i++)
152                s[i].insert(i);
153        p   = f.RevSparseJac(m, s);
154
155        // Use forward mode to compute columns 0 and 2
156        // (make sure order of rows and columns does not matter)
157        CPPAD_TESTVECTOR(size_t) r(3), c(3);
158        VectorBase jac(3);
159        r[0] = 2; c[0] = 0;
160        r[1] = 1; c[1] = 2;
161        r[2] = 0; c[2] = 0;
162        CppAD::sparse_jacobian_work work;
163        size_t n_sweep = f.SparseJacobianForward(x, p, r, c, jac, work);
164        for(k = 0; k < 3; k++)
165        {       ell = r[k] * n + c[k];
166                ok &=  NearEqual(check[ell], jac[k], eps, eps);
167        }
168        ok &= (n_sweep == 1);
169
170        // Use reverse mode to compute rows 0 and 1
171        // (make sure order of rows and columns does not matter)
172        r.resize(4), c.resize(4); jac.resize(4);
173        r[0] = 0; c[0] = 0;
174        r[1] = 1; c[1] = 2;
175        r[2] = 0; c[2] = 1;
176        r[3] = 1; c[3] = 3;
177        work.clear();
178        n_sweep = f.SparseJacobianReverse(x, p, r, c, jac, work);
179        for(k = 0; k < 4; k++)
180        {       ell = r[k] * n + c[k];
181                ok &=  NearEqual(check[ell], jac[k], eps, eps);
182        }
183        ok &= (n_sweep == 1);
184
185        return ok;
186}
187template <class VectorBase, class VectorBool> 
188bool rc_bool(void)
189{       bool ok = true;
190        using CppAD::AD;
191        using CppAD::NearEqual;
192        size_t j, k, ell;
193        double eps = 10. * CppAD::epsilon<double>();
194
195        // domain space vector
196        size_t n = 4;
197        CPPAD_TESTVECTOR(AD<double>)  X(n);
198        for(j = 0; j < n; j++)
199                X[j] = AD<double> (0);
200
201        // declare independent variables and starting recording
202        CppAD::Independent(X);
203
204        size_t m = 3;
205        CPPAD_TESTVECTOR(AD<double>)  Y(m);
206        Y[0] = 1.0*X[0] + 2.0*X[1];
207        Y[1] = 3.0*X[2] + 4.0*X[3];
208        Y[2] = 5.0*X[0] + 6.0*X[1] + 7.0*X[3]*X[3]/2.;
209
210        // create f: X -> Y and stop tape recording
211        CppAD::ADFun<double> f(X, Y);
212
213        // new value for the independent variable vector
214        VectorBase x(n);
215        for(j = 0; j < n; j++)
216                x[j] = double(j);
217
218        // Jacobian of y
219        /*
220              [ 1 2 0 0    ]
221        jac = [ 0 0 3 4    ]
222              [ 5 6 0 7*x_3]
223        */
224        VectorBase check(m * n);
225        check[0] = 1.; check[1] = 2.; check[2]  = 0.; check[3]  = 0.;
226        check[4] = 0.; check[5] = 0.; check[6]  = 3.; check[7]  = 4.;
227        check[8] = 5.; check[9] = 6.; check[10] = 0.; check[11] = 7.*x[3];
228        VectorBool s(m * n);
229        s[0] = true;   s[1] = true;   s[2] = false;   s[3] = false;
230        s[4] = false;  s[5] = false;  s[6] = true;    s[7] = true;
231        s[8] = true;   s[9] = true;  s[10] = false;  s[11] = true;
232
233        // Use forward mode to compute columns 0 and 2
234        // (make sure order of rows and columns does not matter)
235        CPPAD_TESTVECTOR(size_t) r(3), c(3);
236        VectorBase jac(3);
237        r[0] = 2; c[0] = 0;
238        r[1] = 1; c[1] = 2;
239        r[2] = 0; c[2] = 0;
240        CppAD::sparse_jacobian_work work;
241        size_t n_sweep = f.SparseJacobianForward(x, s, r, c, jac, work);
242        for(k = 0; k < 3; k++)
243        {       ell = r[k] * n + c[k];
244                ok &=  NearEqual(check[ell], jac[k], eps, eps);
245        }
246        ok &= (n_sweep == 1);
247
248        // Use reverse mode to compute rows 0 and 1
249        // (make sure order of rows and columns does not matter)
250        r.resize(4), c.resize(4); jac.resize(4);
251        r[0] = 0; c[0] = 0;
252        r[1] = 1; c[1] = 2;
253        r[2] = 0; c[2] = 1;
254        r[3] = 1; c[3] = 3;
255        work.clear();
256        n_sweep = f.SparseJacobianReverse(x, s, r, c, jac, work);
257        for(k = 0; k < 4; k++)
258        {       ell = r[k] * n + c[k];
259                ok &=  NearEqual(check[ell], jac[k], eps, eps);
260        }
261        ok &= (n_sweep == 1);
262
263        return ok;
264}
265
266
267template <class VectorBase, class VectorBool> 
268bool reverse_bool(void)
269{       bool ok = true;
270        using CppAD::AD;
271        using CppAD::NearEqual;
272        size_t i, j, k;
273
274        // domain space vector
275        size_t n = 4;
276        CPPAD_TESTVECTOR(AD<double>)  X(n);
277        for(j = 0; j < n; j++)
278                X[j] = AD<double> (0);
279
280        // declare independent variables and starting recording
281        CppAD::Independent(X);
282
283        size_t m = 3;
284        CPPAD_TESTVECTOR(AD<double>)  Y(m);
285        Y[0] = 1.0*X[0] + 2.0*X[1];
286        Y[1] = 3.0*X[2] + 4.0*X[3];
287        Y[2] = 5.0*X[0] + 6.0*X[1] + 7.0*X[2] + 8.0*X[3]*X[3]/2.;
288
289        // create f: X -> Y and stop tape recording
290        CppAD::ADFun<double> f(X, Y);
291
292        // new value for the independent variable vector
293        VectorBase x(n);
294        for(j = 0; j < n; j++)
295                x[j] = double(j);
296
297        // Jacobian of y without sparsity pattern
298        VectorBase jac(m * n);
299        jac = f.SparseJacobian(x);
300        /*
301              [ 1 2 0 0    ]
302        jac = [ 0 0 3 4    ]
303              [ 5 6 7 8*x_3]
304        */
305        VectorBase check(m * n);
306        check[0] = 1.; check[1] = 2.; check[2]  = 0.; check[3]  = 0.;
307        check[4] = 0.; check[5] = 0.; check[6]  = 3.; check[7]  = 4.;
308        check[8] = 5.; check[9] = 6.; check[10] = 7.; check[11] = 8.*x[3];
309        for(k = 0; k < 12; k++)
310                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
311
312        // test passing sparsity pattern
313        VectorBool s(m * m);
314        VectorBool p(m * n);
315        for(i = 0; i < m; i++)
316        {       for(k = 0; k < m; k++)
317                        s[i * m + k] = false;
318                s[i * m + i] = true;
319        }
320        p   = f.RevSparseJac(m, s);
321        jac = f.SparseJacobian(x);
322        for(k = 0; k < 12; k++)
323                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
324
325        return ok;
326}
327
328template <class VectorBase, class VectorSet> 
329bool reverse_set(void)
330{       bool ok = true;
331        using CppAD::AD;
332        using CppAD::NearEqual;
333        size_t i, j, k;
334
335        // domain space vector
336        size_t n = 4;
337        CPPAD_TESTVECTOR(AD<double>)  X(n);
338        for(j = 0; j < n; j++)
339                X[j] = AD<double> (0);
340
341        // declare independent variables and starting recording
342        CppAD::Independent(X);
343
344        size_t m = 3;
345        CPPAD_TESTVECTOR(AD<double>)  Y(m);
346        Y[0] = X[0] + X[1];
347        Y[1] = X[2] + X[3];
348        Y[2] = X[0] + X[1] + X[2] + X[3] * X[3] / 2.;
349
350        // create f: X -> Y and stop tape recording
351        CppAD::ADFun<double> f(X, Y);
352
353        // new value for the independent variable vector
354        VectorBase x(n);
355        for(j = 0; j < n; j++)
356                x[j] = double(j);
357
358        // Jacobian of y without sparsity pattern
359        VectorBase jac(m * n);
360        jac = f.SparseJacobian(x);
361        /*
362              [ 1 1 0 0  ]
363        jac = [ 0 0 1 1  ]
364              [ 1 1 1 x_3]
365        */
366        VectorBase check(m * n);
367        check[0] = 1.; check[1] = 1.; check[2]  = 0.; check[3]  = 0.;
368        check[4] = 0.; check[5] = 0.; check[6]  = 1.; check[7]  = 1.;
369        check[8] = 1.; check[9] = 1.; check[10] = 1.; check[11] = x[3];
370        for(k = 0; k < 12; k++)
371                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
372
373        // test passing sparsity pattern
374        VectorSet s(m), p(m);
375        for(i = 0; i < m; i++)
376                s[i].insert(i);
377        p   = f.RevSparseJac(m, s);
378        jac = f.SparseJacobian(x);
379        for(k = 0; k < 12; k++)
380                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
381
382        return ok;
383}
384
385template <class VectorBase, class VectorBool> 
386bool forward_bool(void)
387{       bool ok = true;
388        using CppAD::AD;
389        using CppAD::NearEqual;
390        size_t j, k;
391
392        // domain space vector
393        size_t n = 3;
394        CPPAD_TESTVECTOR(AD<double>)  X(n);
395        for(j = 0; j < n; j++)
396                X[j] = AD<double> (0);
397
398        // declare independent variables and starting recording
399        CppAD::Independent(X);
400
401        size_t m = 4;
402        CPPAD_TESTVECTOR(AD<double>)  Y(m);
403        Y[0] = X[0] + X[2];
404        Y[1] = X[0] + X[2];
405        Y[2] = X[1] + X[2];
406        Y[3] = X[1] + X[2] * X[2] / 2.;
407
408        // create f: X -> Y and stop tape recording
409        CppAD::ADFun<double> f(X, Y);
410
411        // new value for the independent variable vector
412        VectorBase x(n);
413        for(j = 0; j < n; j++)
414                x[j] = double(j);
415
416        // Jacobian of y without sparsity pattern
417        VectorBase jac(m * n);
418        jac = f.SparseJacobian(x);
419        /*
420              [ 1 0 1   ]
421        jac = [ 1 0 1   ]
422              [ 0 1 1   ]
423              [ 0 1 x_2 ]
424        */
425        VectorBase check(m * n);
426        check[0] = 1.; check[1]  = 0.; check[2]  = 1.; 
427        check[3] = 1.; check[4]  = 0.; check[5]  = 1.;
428        check[6] = 0.; check[7]  = 1.; check[8]  = 1.; 
429        check[9] = 0.; check[10] = 1.; check[11] = x[2];
430        for(k = 0; k < 12; k++)
431                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
432
433        // test passing sparsity pattern
434        VectorBool r(n * n);
435        VectorBool p(m * n);
436        for(j = 0; j < n; j++)
437        {       for(k = 0; k < n; k++)
438                        r[j * n + k] = false;
439                r[j * n + j] = true;
440        }
441        p   = f.ForSparseJac(n, r);
442        jac = f.SparseJacobian(x);
443        for(k = 0; k < 12; k++)
444                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
445
446        return ok;
447}
448
449template <class VectorBase, class VectorSet> 
450bool forward_set(void)
451{       bool ok = true;
452        using CppAD::AD;
453        using CppAD::NearEqual;
454        size_t j, k;
455
456        // domain space vector
457        size_t n = 3;
458        CPPAD_TESTVECTOR(AD<double>)  X(n);
459        for(j = 0; j < n; j++)
460                X[j] = AD<double> (0);
461
462        // declare independent variables and starting recording
463        CppAD::Independent(X);
464
465        size_t m = 4;
466        CPPAD_TESTVECTOR(AD<double>)  Y(m);
467        Y[0] = X[0] + X[2];
468        Y[1] = X[0] + X[2];
469        Y[2] = X[1] + X[2];
470        Y[3] = X[1] + X[2] * X[2] / 2.;
471
472        // create f: X -> Y and stop tape recording
473        CppAD::ADFun<double> f(X, Y);
474
475        // new value for the independent variable vector
476        VectorBase x(n);
477        for(j = 0; j < n; j++)
478                x[j] = double(j);
479
480        // Jacobian of y without sparsity pattern
481        VectorBase jac(m * n);
482        jac = f.SparseJacobian(x);
483        /*
484              [ 1 0 1   ]
485        jac = [ 1 0 1   ]
486              [ 0 1 1   ]
487              [ 0 1 x_2 ]
488        */
489        VectorBase check(m * n);
490        check[0] = 1.; check[1]  = 0.; check[2]  = 1.; 
491        check[3] = 1.; check[4]  = 0.; check[5]  = 1.;
492        check[6] = 0.; check[7]  = 1.; check[8]  = 1.; 
493        check[9] = 0.; check[10] = 1.; check[11] = x[2];
494        for(k = 0; k < 12; k++)
495                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
496
497        // test passing sparsity pattern
498        VectorSet r(n), p(m);
499        for(j = 0; j < n; j++)
500                r[j].insert(j);
501        p   = f.ForSparseJac(n, r);
502        jac = f.SparseJacobian(x);
503        for(k = 0; k < 12; k++)
504                ok &=  NearEqual(check[k], jac[k], 1e-10, 1e-10 );
505
506        return ok;
507}
508
509bool multiple_of_n_bit(void)
510{       bool ok = true;
511        using CppAD::AD;
512        using CppAD::vector;
513        size_t i, j;
514
515        // should be the same as the corresponding typedef in
516        // cppad/local/sparse_pack.hpp
517        typedef size_t Pack;
518
519        // number of bits per packed value
520        size_t n_bit = std::numeric_limits<Pack>::digits;
521
522        // check case where number of variables is equal to n_bit
523        vector< AD<double> > x(n_bit);
524        vector< AD<double> > y(n_bit);
525
526        // create an AD function with domain and range dimension equal to n_bit
527        CppAD::Independent(x);
528        for(i = 0; i < n_bit; i++)
529                y[i] = x[n_bit - i - 1];
530        CppAD::ADFun<double> f(x, y);
531
532        // Jacobian sparsity patterns
533        vector<bool> r(n_bit * n_bit);
534        vector<bool> s(n_bit * n_bit);
535        for(i = 0; i < n_bit; i++)
536        {       for(j = 0; j < n_bit; j++)
537                        r[ i * n_bit + j ] = (i == j);
538        }
539        s = f.ForSparseJac(n_bit, r);
540
541        // check the result
542        for(i = 0; i < n_bit; i++)
543        {       for(j = 0; j < n_bit; j++)
544                {       if( i == n_bit - j - 1 )
545                                ok = ok & s[ i * n_bit + j ];
546                        else    ok = ok & (! s[i * n_bit + j] );
547                }
548        }
549
550        return ok;
551}
552} // End empty namespace
553# include <vector>
554# include <valarray>
555bool sparse_jacobian(void)
556{       bool ok = true;
557        ok &= rc_tridiagonal();
558        ok &= multiple_of_n_bit();
559        // ---------------------------------------------------------------
560        // vector of bool cases
561        ok &=      rc_bool< CppAD::vector<double>, CppAD::vectorBool   >();
562        ok &= forward_bool< CppAD::vector<double>, CppAD::vector<bool> >();
563        //
564        ok &= reverse_bool< std::vector<double>,   std::vector<bool>   >();
565        ok &=      rc_bool< std::vector<double>,   std::valarray<bool> >();
566        //
567        ok &= forward_bool< std::valarray<double>, CppAD::vectorBool   >();
568        ok &= reverse_bool< std::valarray<double>, CppAD::vector<bool> >();
569        // ---------------------------------------------------------------
570        // vector of set cases
571        typedef std::vector< std::set<size_t> >   std_vector_set;
572        typedef std::valarray< std::set<size_t> > std_valarray_set;
573        typedef CppAD::vector< std::set<size_t> > cppad_vector_set;
574        //
575        ok &=      rc_set< CppAD::vector<double>, std_vector_set   >();
576        ok &= forward_set< std::valarray<double>, std_vector_set   >();
577        //
578        ok &= reverse_set< std::vector<double>,   cppad_vector_set >();
579        ok &=      rc_set< CppAD::vector<double>, cppad_vector_set >();
580        //
581        // According to section 26.3.2.3 of the 1998 C++ standard
582        // a const valarray does not return references to its elements.
583        // ok &= forward_set< std::valarray<double>, std_valarray_set >();
584        // ok &= reverse_set< std::valarray<double>, std_valarray_set >();
585        // ---------------------------------------------------------------
586        //
587        return ok;
588}
Note: See TracBrowser for help on using the repository browser.