source: trunk/test_more/old_usead_1.cpp @ 3740

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

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

commit cd48f5689e158e2758a3bfb57f2e6285b62b8ec1
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 18:54:45 2015 -0700

Run reduce_index on files that were changed today.


reduce_index.py: add some more words to ignore.

commit cf4a8e3a59f06d919ef9a2d911d7c1e462896b39
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 18:35:02 2015 -0700

Change depreciated to deprecated (typo).

commit 37e09901ab78be253fa50004cdc041b042015eeb
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 18:17:29 2015 -0700

Remove invisible white space.

commit 5ffd5cb372caac6c371af5718b8bfc437c8d7f97
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 18:17:09 2015 -0700

Add the deprecation date to some features that were missing it.

commit de437774c18eda2c619cfd8676985113207eb736
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 15:35:52 2015 -0700

whats_new_15.omh: update from previous commit.
test_more.cpp: alphabetize and white space.

commit e8f32f25fe12edf8a8ffa29eb8d2848fef007bd1
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 15:21:17 2015 -0700

  1. Alphabatize list.
  2. remove old_mat_mul.hpp (deprecated).
  3. Correct $cref -> $rref.

commit f205c22e6c4bcabd516214d40e18c4bc4d61a9d1
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 14:44:00 2015 -0700

Remove invisible white space.

commit 1a73e41a8880f8fe351ba4d620cd25055b2a8e8d
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 14:43:43 2015 -0700

Move some deprecated examples from example directory to test_more directory.

commit 2ea5990965e3f8fd00a6393a60b0a59152b41e15
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 14:36:34 2015 -0700

check_user_def.sh: must not include check_user_def.sh in files searched.

commit 26ae3bc6ec918032701b7c1aa37c5b8f9acc6d55
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 13:16:26 2015 -0700

Remove invisible white space.

commit 6cbc289840bc232aff1f676d9ab3b795028aaa36
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 13:16:22 2015 -0700

Add some of the missing dates when certain features were deprecated.

commit 8e5544f467fec5d60eeb1a3739fc024bd5305014
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 09:23:13 2015 -0700

Automatic check of user API preprocessor symbols and its indicated corrections.

File size: 10.3 KB
Line 
1// $Id$
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/*
14$begin old_usead_1.cpp$$
15$spell
16        checkpoint
17        var
18$$
19
20$section Using AD to Compute Atomic Function Derivatives$$
21$mindex inside user checkpoint$$
22
23$head Deprecated 2013-05-27$$
24This example has been deprecated because it is easier to use the
25$cref checkpoint$$ class instead.
26
27$head Purpose$$
28Consider the case where an inner function is used repeatedly in the
29definition of an outer function.
30In this case, it may reduce the number of variables
31$cref/size_var/seq_property/size_var/$$,
32and hence the required memory.
33
34$head Simple Case$$
35This example is the same as $cref old_reciprocal.cpp$$, except that it
36uses AD to compute the
37derivatives needed by an atomic function.
38This is a simple example of an inner function, and hence not really
39useful for the purpose above;
40see $cref old_usead_2.cpp$$ for a more complete example.
41
42$code
43$verbatim%test_more/old_usead_1.cpp%0%// BEGIN C++%// END C++%1%$$
44$$
45
46$end
47*/
48// BEGIN C++
49# include <cppad/cppad.hpp>
50
51namespace { // Begin empty namespace
52        using CppAD::AD;
53        using CppAD::ADFun;
54        using CppAD::vector;
55
56        // ----------------------------------------------------------------------
57        // function that computes reciprocal
58        ADFun<double>* r_ptr_;
59        void create_r(void)
60        {       vector< AD<double> > ax(1), ay(1);
61                ax[0]  = 1;
62                CppAD::Independent(ax);
63                ay[0]  = 1.0 / ax[0];
64                r_ptr_ = new ADFun<double>(ax, ay);
65        }
66        void destroy_r(void)
67        {       delete r_ptr_;
68                r_ptr_ = CPPAD_NULL;
69        }
70
71        // ----------------------------------------------------------------------
72        // forward mode routine called by CppAD
73        bool reciprocal_forward(
74                size_t                   id ,
75                size_t                    k ,
76                size_t                    n ,
77                size_t                    m ,
78                const vector<bool>&      vx ,
79                vector<bool>&            vy ,
80                const vector<double>&    tx ,
81                vector<double>&          ty
82        )
83        {       assert( id == 0 );
84                assert( n == 1 );
85                assert( m == 1 );
86                assert( k == 0 || vx.size() == 0 );
87                bool ok = true;
88                vector<double> x_q(1), y_q(1);
89
90                // check for special case
91                if( vx.size() > 0 )
92                        vy[0] = vx[0];
93
94                // make sure r_ has proper lower order Taylor coefficients stored
95                // then compute ty[k]
96                for(size_t q = 0; q <= k; q++)
97                {       x_q[0] = tx[q];
98                        y_q    = r_ptr_->Forward(q, x_q);
99                        if( q == k )
100                                ty[k] = y_q[0];
101                        assert( q == k || ty[q] == y_q[0] );
102                }
103                return ok;
104        }
105        // ----------------------------------------------------------------------
106        // reverse mode routine called by CppAD
107        bool reciprocal_reverse(
108                size_t                   id ,
109                size_t                    k ,
110                size_t                    n ,
111                size_t                    m ,
112                const vector<double>&    tx ,
113                const vector<double>&    ty ,
114                vector<double>&          px ,
115                const vector<double>&    py
116        )
117        {       assert( id == 0 );
118                assert( n == 1 );
119                assert( m == 1 );
120                bool ok = true;
121                vector<double> x_q(1), w(k+1), dw(k+1);
122
123                // make sure r_ has proper forward mode coefficients
124                size_t q;
125                for(q = 0; q <= k; q++)
126                {       x_q[0] = tx[q];
127# ifdef NDEBUG
128                        r_ptr_->Forward(q, x_q);
129# else
130                        vector<double> y_q(1);
131                        y_q    = r_ptr_->Forward(q, x_q);
132                        assert( ty[q] == y_q[0] );
133# endif
134                }
135                for(q = 0; q <=k; q++)
136                        w[q] = py[q];
137                dw = r_ptr_->Reverse(k+1, w);
138                for(q = 0; q <=k; q++)
139                        px[q] = dw[q];
140
141                return ok;
142        }
143        // ----------------------------------------------------------------------
144        // forward Jacobian sparsity routine called by CppAD
145        bool reciprocal_for_jac_sparse(
146                size_t                               id ,
147                size_t                                n ,
148                size_t                                m ,
149                size_t                                p ,
150                const vector< std::set<size_t> >&     r ,
151                vector< std::set<size_t> >&           s )
152        {       assert( id == 0 );
153                assert( n == 1 );
154                assert( m == 1 );
155                bool ok = true;
156
157                vector< std::set<size_t> > R(1), S(1);
158                R[0] = r[0];
159                S = r_ptr_->ForSparseJac(p, R);
160                s[0] = S[0];
161
162                return ok;
163        }
164        // ----------------------------------------------------------------------
165        // reverse Jacobian sparsity routine called by CppAD
166        bool reciprocal_rev_jac_sparse(
167                size_t                               id ,
168                size_t                                n ,
169                size_t                                m ,
170                size_t                                p ,
171                vector< std::set<size_t> >&           r ,
172                const vector< std::set<size_t> >&     s )
173        {
174                assert( id == 0 );
175                assert( n == 1 );
176                assert( m == 1 );
177                bool ok = true;
178
179                vector< std::set<size_t> > R(p), S(p);
180                size_t q;
181                for(q = 0; q < p; q++)
182                        S[q] = s[q];
183                R = r_ptr_->RevSparseJac(p, S);
184                for(q = 0; q < p; q++)
185                        r[q] = R[q];
186
187                return ok;
188        }
189        // ----------------------------------------------------------------------
190        // reverse Hessian sparsity routine called by CppAD
191        bool reciprocal_rev_hes_sparse(
192                size_t                               id ,
193                size_t                                n ,
194                size_t                                m ,
195                size_t                                p ,
196                const vector< std::set<size_t> >&     r ,
197                const vector<bool>&                   s ,
198                vector<bool>&                         t ,
199                const vector< std::set<size_t> >&     u ,
200                vector< std::set<size_t> >&           v )
201        {       // Can just return false if not use RevSparseHes.
202                assert( id == 0 );
203                assert( n == 1 );
204                assert( m == 1 );
205                bool ok = true;
206
207                // compute sparsity pattern for T(x) = S(x) * f'(x)
208                vector<bool> T(1), S(1);
209                S[0]   = s[0];
210                T      = r_ptr_->RevSparseJac(1, S);
211                t[0]   = T[0];
212
213                // compute sparsity pattern for A(x) = U(x)^T * f'(x)
214                vector<bool> Ut(p), A(p);
215                size_t q;
216                for(q = 0; q < p; q++)
217                        Ut[q] = false;
218                std::set<size_t>::iterator itr;
219                for(itr = u[0].begin(); itr != u[0].end(); itr++)
220                        Ut[*itr] = true;
221                A = r_ptr_-> RevSparseJac(p, Ut);
222
223                // compute sparsity pattern for H(x) = R^T * (S * F)''(x)
224                vector<bool> H(p), R(n);
225                for(q = 0; q < p; q++)
226                        R[q] = false;
227                for(itr = r[0].begin(); itr != r[0].end(); itr++)
228                        R[*itr] = true;
229                r_ptr_->ForSparseJac(p, R);
230                H = r_ptr_->RevSparseHes(p, S);
231
232                // compute sparsity pattern for V(x) = A(x)^T + H(x)^T
233                v[0].clear();
234                for(q = 0; q < p; q++)
235                        if( A[q] | H[q] )
236                                v[0].insert(q);
237
238                return ok;
239        }
240        // ---------------------------------------------------------------------
241        // Declare the AD<double> routine reciprocal(id, ax, ay)
242        CPPAD_USER_ATOMIC(
243                reciprocal                 ,
244                CppAD::vector              ,
245                double                     ,
246                reciprocal_forward         ,
247                reciprocal_reverse         ,
248                reciprocal_for_jac_sparse  ,
249                reciprocal_rev_jac_sparse  ,
250                reciprocal_rev_hes_sparse
251        )
252} // End empty namespace
253
254bool old_usead_1(void)
255{       bool ok = true;
256        using CppAD::NearEqual;
257        double eps = 10. * CppAD::numeric_limits<double>::epsilon();
258
259        // --------------------------------------------------------------------
260        // Create the ADFun<doulbe> r_
261        create_r();
262
263        // --------------------------------------------------------------------
264        // Create the function f(x)
265        //
266        // domain space vector
267        size_t n  = 1;
268        double  x0 = 0.5;
269        vector< AD<double> > ax(n);
270        ax[0]     = x0;
271
272        // declare independent variables and start tape recording
273        CppAD::Independent(ax);
274
275        // range space vector
276        size_t m = 1;
277        vector< AD<double> > ay(m);
278
279        // call user function and store reciprocal(x) in au[0]
280        vector< AD<double> > au(m);
281        size_t id = 0;           // not used
282        reciprocal(id, ax, au); // u = 1 / x
283
284        // call user function and store reciprocal(u) in ay[0]
285        reciprocal(id, au, ay); // y = 1 / u = x
286
287        // create f: x -> y and stop tape recording
288        ADFun<double> f;
289        f.Dependent(ax, ay);  // f(x) = x
290
291        // --------------------------------------------------------------------
292        // Check function value results
293        //
294        // check function value
295        double check = x0;
296        ok &= NearEqual( Value(ay[0]) , check,  eps, eps);
297
298        // check zero order forward mode
299        size_t q;
300        vector<double> x_q(n), y_q(m);
301        q      = 0;
302        x_q[0] = x0;
303        y_q    = f.Forward(q, x_q);
304        ok &= NearEqual(y_q[0] , check,  eps, eps);
305
306        // check first order forward mode
307        q      = 1;
308        x_q[0] = 1;
309        y_q    = f.Forward(q, x_q);
310        check  = 1.;
311        ok &= NearEqual(y_q[0] , check,  eps, eps);
312
313        // check second order forward mode
314        q      = 2;
315        x_q[0] = 0;
316        y_q    = f.Forward(q, x_q);
317        check  = 0.;
318        ok &= NearEqual(y_q[0] , check,  eps, eps);
319
320        // --------------------------------------------------------------------
321        // Check reverse mode results
322        //
323        // third order reverse mode
324        q     = 3;
325        vector<double> w(m), dw(n * q);
326        w[0]  = 1.;
327        dw    = f.Reverse(q, w);
328        check = 1.;
329        ok &= NearEqual(dw[0] , check,  eps, eps);
330        check = 0.;
331        ok &= NearEqual(dw[1] , check,  eps, eps);
332        ok &= NearEqual(dw[2] , check,  eps, eps);
333
334        // --------------------------------------------------------------------
335        // forward mode sparstiy pattern
336        size_t p = n;
337        CppAD::vectorBool r1(n * p), s1(m * p);
338        r1[0] = true;          // compute sparsity pattern for x[0]
339        s1    = f.ForSparseJac(p, r1);
340        ok  &= s1[0] == true;  // f[0] depends on x[0]
341
342        // --------------------------------------------------------------------
343        // reverse mode sparstiy pattern
344        q = m;
345        CppAD::vectorBool s2(q * m), r2(q * n);
346        s2[0] = true;          // compute sparsity pattern for f[0]
347        r2    = f.RevSparseJac(q, s2);
348        ok  &= r2[0] == true;  // f[0] depends on x[0]
349
350        // --------------------------------------------------------------------
351        // Hessian sparsity (using previous ForSparseJac call)
352        CppAD::vectorBool s3(m), h(p * n);
353        s3[0] = true;        // compute sparsity pattern for f[0]
354        h     = f.RevSparseJac(p, s3);
355        ok  &= h[0] == true; // second partial of f[0] w.r.t. x[0] may be non-zero
356
357        // -----------------------------------------------------------------
358        // Free all memory associated with the object r_ptr
359        destroy_r();
360
361        // -----------------------------------------------------------------
362        // Free all temporary work space associated with old_atomic objects.
363        // (If there are future calls to user atomic functions, they will
364        // create new temporary work space.)
365        CppAD::user_atomic<double>::clear();
366
367        return ok;
368}
369// END C++
Note: See TracBrowser for help on using the repository browser.