source: trunk/test_more/checkpoint.cpp @ 3008

Last change on this file since 3008 was 2859, checked in by bradbell, 7 years ago

merge in changes from branches/atomic; see bin/svn_merge.sh

  • Property svn:keywords set to Id
File size: 3.8 KB
Line 
1/* $Id: checkpoint.cpp 2859 2013-05-28 06:03:21Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 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
40bool checkpoint(void)
41{       bool ok = true;
42        using CppAD::checkpoint;
43        using CppAD::ADFun;
44        using CppAD::NearEqual;
45        size_t i, j, k, n = 4, ell = n-1 , m = ell + 1;
46        double eps = 10. * std::numeric_limits<double>::epsilon();
47
48        // checkpoint version of the function F(x)
49        ADVector ax(n), ay(ell), az(m);
50        for(j = 0; j < n; j++)
51                ax[j] = double(j);
52        checkpoint<double> f_check("f_check", f_algo, ax, ay); 
53        checkpoint<double> g_check("g_check", g_algo, ay, az); 
54
55        // Record a version of z = g[f(x)] without checkpointing
56        Independent(ax);
57        f_algo(ax, ay);
58        g_algo(ay, az);
59        ADFun<double> check_not(ax, az);
60
61        // Record a version of z = g[f(x)] with checkpointing
62        Independent(ax);
63        f_check(ax, ay);
64        g_check(ay, az);
65        ADFun<double> check_yes(ax, az);
66
67        // compare forward mode results for orders 0, 1, 2
68        size_t p = 2;
69        CPPAD_TESTVECTOR(double) x_p(n*(p+1)), z_not(m*(p+1)), z_yes(m*(p+1));
70        for(j = 0; j < n; j++)
71        {       for(k = 0; k <= p; k++)
72                        x_p[ j * (p+1) + k ] = 1.0 / (p + 1 - k);
73        }
74        z_not = check_not.Forward(p, x_p);
75        z_yes = check_yes.Forward(p, x_p);
76        for(i = 0; i < m; i++)
77        {       for(k = 0; k <= p; k++)
78                {       double zik_not = z_not[ i * (p+1) + k];
79                        double zik_yes = z_yes[ i * (p+1) + k];
80                        ok &= NearEqual(zik_not, zik_yes, eps, eps);
81                }
82        }
83
84        // compare reverse mode results
85        CPPAD_TESTVECTOR(double) w(m*(p+1)), dw_not(n*(p+1)), dw_yes(n*(p+1));
86        dw_not = check_not.Reverse(p+1, w);
87        dw_yes = check_yes.Reverse(p+1, w);
88        for(j = 0; j < n; j++)
89        {       for(k = 0; k <= p; k++)
90                {       double dwjk_not = dw_not[ j * (p+1) + k];
91                        double dwjk_yes = dw_yes[ j * (p+1) + k];
92                        ok &= NearEqual(dwjk_not, dwjk_yes, eps, eps);
93                }
94        }
95
96        // mix sparsity so test both cases
97        f_check.option( CppAD::atomic_base<double>::bool_sparsity_enum );
98        g_check.option( CppAD::atomic_base<double>::set_sparsity_enum );
99
100        // compare forward mode Jacobian sparsity patterns
101        size_t q = n - 1;
102        CppAD::vector< std::set<size_t> > r(n), s_not(m), s_yes(m);
103        for(j = 0; j < n; j++)
104        {       if( j < q )
105                        r[j].insert(j);
106                else
107                {       r[j].insert(0);
108                        r[j].insert(1);
109                }
110        }
111        s_not = check_not.ForSparseJac(q, r);
112        s_yes = check_yes.ForSparseJac(q, r);
113        for(i = 0; i < m; i++)
114                ok &= s_not[i] == s_yes[i];
115
116        // compare reverse mode Jacobian sparsity patterns
117        CppAD::vector< std::set<size_t> > s(m), r_not(m), r_yes(m);
118        for(i = 0; i < m; i++)
119                s[i].insert(i);
120        r_not = check_not.RevSparseJac(m, s);
121        r_yes = check_yes.RevSparseJac(m, s);
122        for(i = 0; i < m; i++)
123                ok &= s_not[i] == s_yes[i];
124       
125
126        // compare reverse mode Hessian sparsity patterns
127        CppAD::vector< std::set<size_t> > s_one(1), h_not(q), h_yes(q); 
128        for(i = 0; i < m; i++)
129                s_one[0].insert(i);
130        h_not = check_not.RevSparseHes(q, s_one);
131        h_yes = check_yes.RevSparseHes(q, s_one);
132        for(i = 0; i < q; i++)
133                ok &= h_not[i] == h_yes[i];
134       
135        checkpoint<double>::clear();
136        return ok;
137}
138// END C++
Note: See TracBrowser for help on using the repository browser.