source: trunk/test_more/checkpoint.cpp @ 3793

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

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: 08d0cfc73520f8d330289da641b81d557128d048
end hash code: c4559d5e01e1b0f09943490dd84449557eced25d

commit c4559d5e01e1b0f09943490dd84449557eced25d
Author: Brad Bell <bradbell@…>
Date: Thu Jan 21 15:27:44 2016 -0700

checkpoint.cpp: fix initialization error in test (found by valgrind).

  • Property svn:keywords set to Id
File size: 5.1 KB
Line 
1// $Id: checkpoint.cpp 3784 2016-01-22 10:46:29Z bradbell $
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-16 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                for(i = 0; i < m; i++)
86                {       for(k = 0; k <= p; k++)
87                                w[ i * (p+1) + k ] = 2.0 / (p + 1 - k );
88                }
89                dw_not = check_not.Reverse(p+1, w);
90                dw_yes = check_yes.Reverse(p+1, w);
91                for(j = 0; j < n; j++)
92                {       for(k = 0; k <= p; k++)
93                        {       double dwjk_not = dw_not[ j * (p+1) + k];
94                                double dwjk_yes = dw_yes[ j * (p+1) + k];
95                                ok &= NearEqual(dwjk_not, dwjk_yes, eps, eps);
96                        }
97                }
98
99                // mix sparsity so test both cases
100                f_check.option( CppAD::atomic_base<double>::bool_sparsity_enum );
101                g_check.option( CppAD::atomic_base<double>::set_sparsity_enum );
102
103                // compare forward mode Jacobian sparsity patterns
104                size_t q = n - 1;
105                CppAD::vector< std::set<size_t> > r(n), s_not(m), s_yes(m);
106                for(j = 0; j < n; j++)
107                {       if( j < q )
108                                r[j].insert(j);
109                        else
110                        {       r[j].insert(0);
111                                r[j].insert(1);
112                        }
113                }
114                s_not = check_not.ForSparseJac(q, r);
115                s_yes = check_yes.ForSparseJac(q, r);
116                for(i = 0; i < m; i++)
117                        ok &= s_not[i] == s_yes[i];
118
119                // compare reverse mode Jacobian sparsity patterns
120                CppAD::vector< std::set<size_t> > s(m), r_not(m), r_yes(m);
121                for(i = 0; i < m; i++)
122                        s[i].insert(i);
123                r_not = check_not.RevSparseJac(m, s);
124                r_yes = check_yes.RevSparseJac(m, s);
125                for(i = 0; i < m; i++)
126                        ok &= s_not[i] == s_yes[i];
127
128
129                // compare reverse mode Hessian sparsity patterns
130                CppAD::vector< std::set<size_t> > s_one(1), h_not(q), h_yes(q);
131                for(i = 0; i < m; i++)
132                        s_one[0].insert(i);
133                h_not = check_not.RevSparseHes(q, s_one);
134                h_yes = check_yes.RevSparseHes(q, s_one);
135                for(i = 0; i < q; i++)
136                        ok &= h_not[i] == h_yes[i];
137
138                checkpoint<double>::clear();
139                return ok;
140        }
141
142        bool h_algo(const ADVector& ax, ADVector& ay)
143        {       ay[0] = ax[0];
144                ay[1] = ax[1] + ax[2];
145                return true;
146        }
147        bool test_two(void)
148        {       bool ok = true;
149                using CppAD::checkpoint;
150                using CppAD::ADFun;
151                using CppAD::NearEqual;
152
153                // checkpoint version of H(x)
154                size_t m = 2;
155                size_t n = 3;
156                ADVector ax(n), ay(m);
157                for(size_t j = 0; j < n; j++)
158                        ax[j] = double(j);
159                checkpoint<double> h_check("h_check", h_algo, ax, ay);
160
161                // record function using h_check
162                Independent(ax);
163                h_check(ax, ay);
164                ADFun<double> h(ax, ay);
165
166                for(size_t k = 0; k < 3; k++)
167                {       if( k == 0 )
168                                h_check.option(CppAD::atomic_base<double>::pack_sparsity_enum);
169                        if( k == 1 )
170                                h_check.option(CppAD::atomic_base<double>::bool_sparsity_enum);
171                        if( k == 2 )
172                                h_check.option(CppAD::atomic_base<double>::set_sparsity_enum);
173
174                        // compute sparsity pattern h_1(x) = x[1] + x[2]
175                        CppAD::vector< std::set<size_t> > r(1), s(1);
176                        r[0].insert(1);
177                        s = h.RevSparseJac(1, r);
178
179                        // check result
180                        std::set<size_t> check;
181                        check.insert(1);
182                        check.insert(2);
183                        ok &= s[0] == check;
184                }
185
186                return ok;
187        }
188}
189
190bool checkpoint(void)
191{       bool ok = true;
192        ok  &= test_one();
193        ok  &= test_two();
194        return ok;
195}
196// END C++
Note: See TracBrowser for help on using the repository browser.