source: trunk/test_more/sparse_jacobian.cpp @ 2401

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

Merge in changes from branches/sparse.

  1. Update version number to today.
  2. Add sparse return user API interfaces to SparseJacobian? and SparseHessian?.
  3. Include new interface and imporve sparse Hessian and Jacobian examples.
  4. More testing of sparse Hessian and Jacobians.
  5. Add is_element to sparse_set and sparse_pack.

.: merge information.
search.sh: search cppad/speed directory, drop makefile.in from results.
svn_merge.sh: parameters for this merge.
parallel_ad.hpp: initilize is_element statics.
sparse_evaluate.hpp: fix typo in documentation.
glossary.omh: convert to modern omhelp, improve AD Function entry.
link_sparse_hessian.cpp: fix typo in documentation.

  • Property svn:keywords set to Id
File size: 15.1 KB
Line 
1/* $Id: sparse_jacobian.cpp 2401 2012-05-24 16:30:25Z 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                    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
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_TEST_VECTOR< AD<double> >  X(n);
30        CPPAD_TEST_VECTOR<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_TEST_VECTOR< AD<double> >  Y(m);
39        CPPAD_TEST_VECTOR<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_TEST_VECTOR<size_t> r(K), c(K);
75        CPPAD_TEST_VECTOR<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_TEST_VECTOR< 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_TEST_VECTOR< 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_TEST_VECTOR<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_TEST_VECTOR< 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_TEST_VECTOR< 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_TEST_VECTOR<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_TEST_VECTOR< 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_TEST_VECTOR< 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_TEST_VECTOR< 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_TEST_VECTOR< 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_TEST_VECTOR< 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_TEST_VECTOR< 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_TEST_VECTOR< 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_TEST_VECTOR< 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.