source: trunk/test_more/checkpoint.cpp @ 3713

Last change on this file since 3713 was 3713, checked in by bradbell, 4 years ago

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: fc86286a32cd87ded0475c12f6a8ea5388bc0ccd
end hash code: 630f51c9e1a0037b2302ab026d1f2867b9490e91

commit 630f51c9e1a0037b2302ab026d1f2867b9490e91
Author: Brad Bell <bradbell@…>
Date: Fri Aug 28 05:49:06 2015 -0700

Make checkpoint rev_sparse_jac bool case more like set case.


checkpoint.hpp: fix formating in documentation, set n, m in normal way.
checkpoint.cpp: extend test from set to both bool and set cases.

commit 4f295ecc0efcaead27b701c5e7265ebb6c9e29ba
Author: Brad Bell <bradbell@…>
Date: Thu Aug 27 09:21:25 2015 -0700

automatic freeing of sparsity patterns during forward mode calculations.

  • Property svn:keywords set to Id
File size: 4.9 KB
Line 
1/* $Id: checkpoint.cpp 3713 2015-08-28 13:01:09Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-15 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# include <cppad/cppad.hpp>
14
15namespace {
16        using CppAD::AD;
17        typedef CPPAD_TESTVECTOR(AD<double>) ADVector;
18
19        bool f_algo(const ADVector& x, ADVector& y)
20        {       size_t n = x.size();
21                size_t m = y.size();
22                assert( n == m + 1);
23                for(size_t i = 0; i < m; i++)
24                        y[i] = x[i] * x[i+1];
25                return true;
26        }
27        bool g_algo(const ADVector& y, ADVector& z)
28        {       size_t n = y.size();
29                size_t m = z.size();
30                assert( n + 1 == m );
31                z[0] = 0.0;
32                for(size_t i = 1; i < m; i++)
33                {       z[0] += y[i-1];
34                        z[i]  = y[i-1];
35                }
36                return true;
37        }
38
39        bool test_one(void)
40        {       bool ok = true;
41                using CppAD::checkpoint;
42                using CppAD::ADFun;
43                using CppAD::NearEqual;
44                size_t i, j, k, n = 4, ell = n-1 , m = ell + 1;
45                double eps = 10. * std::numeric_limits<double>::epsilon();
46
47                // checkpoint version of the function F(x)
48                ADVector ax(n), ay(ell), az(m);
49                for(j = 0; j < n; j++)
50                        ax[j] = double(j);
51                checkpoint<double> f_check("f_check", f_algo, ax, ay);
52                checkpoint<double> g_check("g_check", g_algo, ay, az);
53
54                // Record a version of z = g[f(x)] without checkpointing
55                Independent(ax);
56                f_algo(ax, ay);
57                g_algo(ay, az);
58                ADFun<double> check_not(ax, az);
59
60                // Record a version of z = g[f(x)] with checkpointing
61                Independent(ax);
62                f_check(ax, ay);
63                g_check(ay, az);
64                ADFun<double> check_yes(ax, az);
65
66                // compare forward mode results for orders 0, 1, 2
67                size_t p = 2;
68                CPPAD_TESTVECTOR(double) x_p(n*(p+1)), z_not(m*(p+1)), z_yes(m*(p+1));
69                for(j = 0; j < n; j++)
70                {       for(k = 0; k <= p; k++)
71                                x_p[ j * (p+1) + k ] = 1.0 / (p + 1 - k);
72                }
73                z_not = check_not.Forward(p, x_p);
74                z_yes = check_yes.Forward(p, x_p);
75                for(i = 0; i < m; i++)
76                {       for(k = 0; k <= p; k++)
77                        {       double zik_not = z_not[ i * (p+1) + k];
78                                double zik_yes = z_yes[ i * (p+1) + k];
79                                ok &= NearEqual(zik_not, zik_yes, eps, eps);
80                        }
81                }
82
83                // compare reverse mode results
84                CPPAD_TESTVECTOR(double) w(m*(p+1)), dw_not(n*(p+1)), dw_yes(n*(p+1));
85                dw_not = check_not.Reverse(p+1, w);
86                dw_yes = check_yes.Reverse(p+1, w);
87                for(j = 0; j < n; j++)
88                {       for(k = 0; k <= p; k++)
89                        {       double dwjk_not = dw_not[ j * (p+1) + k];
90                                double dwjk_yes = dw_yes[ j * (p+1) + k];
91                                ok &= NearEqual(dwjk_not, dwjk_yes, eps, eps);
92                        }
93                }
94
95                // mix sparsity so test both cases
96                f_check.option( CppAD::atomic_base<double>::bool_sparsity_enum );
97                g_check.option( CppAD::atomic_base<double>::set_sparsity_enum );
98
99                // compare forward mode Jacobian sparsity patterns
100                size_t q = n - 1;
101                CppAD::vector< std::set<size_t> > r(n), s_not(m), s_yes(m);
102                for(j = 0; j < n; j++)
103                {       if( j < q )
104                                r[j].insert(j);
105                        else
106                        {       r[j].insert(0);
107                                r[j].insert(1);
108                        }
109                }
110                s_not = check_not.ForSparseJac(q, r);
111                s_yes = check_yes.ForSparseJac(q, r);
112                for(i = 0; i < m; i++)
113                        ok &= s_not[i] == s_yes[i];
114
115                // compare reverse mode Jacobian sparsity patterns
116                CppAD::vector< std::set<size_t> > s(m), r_not(m), r_yes(m);
117                for(i = 0; i < m; i++)
118                        s[i].insert(i);
119                r_not = check_not.RevSparseJac(m, s);
120                r_yes = check_yes.RevSparseJac(m, s);
121                for(i = 0; i < m; i++)
122                        ok &= s_not[i] == s_yes[i];
123
124
125                // compare reverse mode Hessian sparsity patterns
126                CppAD::vector< std::set<size_t> > s_one(1), h_not(q), h_yes(q);
127                for(i = 0; i < m; i++)
128                        s_one[0].insert(i);
129                h_not = check_not.RevSparseHes(q, s_one);
130                h_yes = check_yes.RevSparseHes(q, s_one);
131                for(i = 0; i < q; i++)
132                        ok &= h_not[i] == h_yes[i];
133
134                checkpoint<double>::clear();
135                return ok;
136        }
137
138        bool h_algo(const ADVector& ax, ADVector& ay)
139        {       ay[0] = ax[0];
140                ay[1] = ax[1] + ax[2];
141                return true;
142        }
143        bool test_two(void)
144        {       bool ok = true;
145                using CppAD::checkpoint;
146                using CppAD::ADFun;
147                using CppAD::NearEqual;
148
149                // checkpoint version of H(x)
150                size_t m = 2;
151                size_t n = 3;
152                ADVector ax(n), ay(m);
153                for(size_t j = 0; j < n; j++)
154                        ax[j] = double(j);
155                checkpoint<double> h_check("h_check", h_algo, ax, ay);
156
157                // record function using h_check
158                Independent(ax);
159                h_check(ax, ay);
160                ADFun<double> h(ax, ay);
161
162                for(size_t k = 0; k < 2; k++)
163                {       if( k == 0 )
164                                h_check.option(CppAD::atomic_base<double>::bool_sparsity_enum);
165                        else
166                                h_check.option(CppAD::atomic_base<double>::set_sparsity_enum);
167
168                        // compute sparsity pattern h_1(x) = x[1] + x[2]
169                        CppAD::vector< std::set<size_t> > r(1), s(1);
170                        r[0].insert(1);
171                        s = h.RevSparseJac(1, r);
172
173                        // check result
174                        ok &= s[0] == std::set<size_t>{1, 2};
175                }
176
177                return ok;
178        }
179}
180
181bool checkpoint(void)
182{       bool ok = true;
183        ok  &= test_one();
184        ok  &= test_two();
185        return ok;
186}
187// END C++
Note: See TracBrowser for help on using the repository browser.