source: trunk/test_more/old_usead_2.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: 14.5 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_2.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
24$head Deprecated 2013-05-27$$
25This example has been deprecated because it is easier to use the
26$cref checkpoint$$ class instead.
27
28$head Purpose$$
29Consider the case where an inner function is used repeatedly in the
30definition of an outer function.
31In this case, it may reduce the number of variables
32$cref/size_var/seq_property/size_var/$$,
33and hence the required memory.
34
35$code
36$verbatim%test_more/old_usead_2.cpp%0%// BEGIN C++%// END C++%1%$$
37$$
38
39$end
40*/
41// BEGIN C++
42# include <cppad/cppad.hpp>
43
44namespace { // Begin empty namespace
45        using CppAD::AD;
46        using CppAD::ADFun;
47        using CppAD::vector;
48
49        // ----------------------------------------------------------------------
50        // ODE for [t, t^2 / 2 ] in form required by Runge45
51        class Fun {
52        public:
53                void Ode(
54                        const AD<double>           &t,
55                        const vector< AD<double> > &z,
56                        vector< AD<double> >       &f)
57                {       assert( z.size() == 2 );
58                        assert( f.size() == 2 );
59                        f[0] =  1.0;
60                        f[1] =  z[0];
61                }
62        };
63
64        // ----------------------------------------------------------------------
65        // Create function that takes on Runge45 step for the ODE above
66        ADFun<double>* r_ptr_;
67        void create_r(void)
68        {       size_t n = 3, m = 2;
69                vector< AD<double> > x(n), zi(m), y(m), e(m);
70                // The value of x does not matter because the operation sequence
71                // does not depend on x.
72                x[0]  = 0.0;  // initial value z_0 (t) at t = ti
73                x[1]  = 0.0;  // initial value z_1 (t) at t = ti
74                x[2]  = 0.1;  // final time for this integration
75                CppAD::Independent(x);
76                zi[0]         = x[0];  // z_0 (t) at t = ti
77                zi[1]         = x[1];  // z_1 (t) at t = ti
78                AD<double> ti = 0.0;   // t does not appear in ODE so does not matter
79                AD<double> tf = x[2];  // final time
80                size_t M      = 3;     // number of Runge45 steps to take
81                Fun F;
82                y             = CppAD::Runge45(F, M, ti, tf, zi, e);
83                r_ptr_        = new ADFun<double>(x, y);
84        }
85        void destroy_r(void)
86        {       delete r_ptr_;
87                r_ptr_ = CPPAD_NULL;
88        }
89
90        // ----------------------------------------------------------------------
91        // forward mode routine called by CppAD
92        bool solve_ode_forward(
93                size_t                   id ,
94                size_t                    k ,
95                size_t                    n ,
96                size_t                    m ,
97                const vector<bool>&      vx ,
98                vector<bool>&            vy ,
99                const vector<double>&    tx ,
100                vector<double>&          ty
101        )
102        {       assert( id == 0 );
103                assert( n == 3 );
104                assert( m == 2 );
105                assert( k == 0 || vx.size() == 0 );
106                bool ok = true;
107                vector<double> xp(n), yp(m);
108                size_t i, j;
109
110                // check for special case
111                if( vx.size() > 0 )
112                {       //Compute r, a Jacobian sparsity pattern.
113                        // Use reverse mode because m < n.
114                        vector< std::set<size_t> > s(m), r(m);
115                        for(i = 0; i < m; i++)
116                                s[i].insert(i);
117                        r = r_ptr_->RevSparseJac(m, s);
118                        std::set<size_t>::const_iterator itr;
119                        for(i = 0; i < m; i++)
120                        {       vy[i] = false;
121                                for(itr = s[i].begin(); itr != s[i].end(); itr++)
122                                {       j = *itr;
123                                        assert( j < n );
124                                        // y[i] depends on the value of x[j]
125                                        // Visual Studio 2013 generates warning without bool below
126                                        vy[i] |= bool( vx[j] );
127                                }
128                        }
129                }
130                // make sure r_ has proper lower order Taylor coefficients stored
131                // then compute ty[k]
132                for(size_t q = 0; q <= k; q++)
133                {       for(j = 0; j < n; j++)
134                                xp[j] = tx[j * (k+1) + q];
135                        yp    = r_ptr_->Forward(q, xp);
136                        if( q == k )
137                        {       for(i = 0; i < m; i++)
138                                        ty[i * (k+1) + q] = yp[i];
139                        }
140# ifndef NDEBUG
141                        else
142                        {       for(i = 0; i < m; i++)
143                                        assert( ty[i * (k+1) + q] == yp[i] );
144                        }
145# endif
146                }
147                // no longer need the Taylor coefficients in r_ptr_
148                // (have to reconstruct them every time)
149                r_ptr_->capacity_order(0);
150                return ok;
151        }
152        // ----------------------------------------------------------------------
153        // reverse mode routine called by CppAD
154        bool solve_ode_reverse(
155                size_t                   id ,
156                size_t                    k ,
157                size_t                    n ,
158                size_t                    m ,
159                const vector<double>&    tx ,
160                const vector<double>&    ty ,
161                vector<double>&          px ,
162                const vector<double>&    py
163        )
164        {       assert( id == 0 );
165                assert( n == 3 );
166                assert( m == 2 );
167                bool ok = true;
168                vector<double> xp(n), w( (k+1) * m ), dw( (k+1) * n );
169
170                // make sure r_ has proper forward mode coefficients
171                size_t i, j, q;
172                for(q = 0; q <= k; q++)
173                {       for(j = 0; j < n; j++)
174                                xp[j] = tx[j * (k+1) + q];
175# ifdef NDEBUG
176                        r_ptr_->Forward(q, xp);
177# else
178                        vector<double> yp(m);
179                        yp = r_ptr_->Forward(q, xp);
180                        for(i = 0; i < m; i++)
181                                assert( ty[i * (k+1) + q] == yp[i] );
182# endif
183                }
184                for(i = 0; i < m; i++)
185                {       for(q = 0; q <=k; q++)
186                                w[ i * (k+1) + q] = py[ i * (k+1) + q];
187                }
188                dw = r_ptr_->Reverse(k+1, w);
189                for(j = 0; j < n; j++)
190                {       for(q = 0; q <=k; q++)
191                                px[ j * (k+1) + q] = dw[ j * (k+1) + q];
192                }
193                // no longer need the Taylor coefficients in r_ptr_
194                // (have to reconstruct them every time)
195                r_ptr_->capacity_order(0);
196
197                return ok;
198        }
199        // ----------------------------------------------------------------------
200        // forward Jacobian sparsity routine called by CppAD
201        bool solve_ode_for_jac_sparse(
202                size_t                               id ,
203                size_t                                n ,
204                size_t                                m ,
205                size_t                                p ,
206                const vector< std::set<size_t> >&     r ,
207                vector< std::set<size_t> >&           s )
208        {       assert( id == 0 );
209                assert( n == 3 );
210                assert( m == 2 );
211                bool ok = true;
212
213                vector< std::set<size_t> > R(n), S(m);
214                for(size_t j = 0; j < n; j++)
215                        R[j] = r[j];
216                S = r_ptr_->ForSparseJac(p, R);
217                for(size_t i = 0; i < m; i++)
218                        s[i] = S[i];
219
220                // no longer need the forward mode sparsity pattern
221                // (have to reconstruct them every time)
222                r_ptr_->size_forward_set(0);
223
224                return ok;
225        }
226        // ----------------------------------------------------------------------
227        // reverse Jacobian sparsity routine called by CppAD
228        bool solve_ode_rev_jac_sparse(
229                size_t                               id ,
230                size_t                                n ,
231                size_t                                m ,
232                size_t                                p ,
233                vector< std::set<size_t> >&           r ,
234                const vector< std::set<size_t> >&     s )
235        {
236                assert( id == 0 );
237                assert( n == 3 );
238                assert( m == 2 );
239                bool ok = true;
240
241                vector< std::set<size_t> > R(p), S(p);
242                std::set<size_t>::const_iterator itr;
243                size_t i;
244                // untranspose s
245                for(i = 0; i < m; i++)
246                {       for(itr = s[i].begin(); itr != s[i].end(); itr++)
247                                S[*itr].insert(i);
248                }
249                R = r_ptr_->RevSparseJac(p, S);
250                // transpose r
251                for(i = 0; i < m; i++)
252                        r[i].clear();
253                for(i = 0; i < p; i++)
254                {       for(itr = R[i].begin(); itr != R[i].end(); itr++)
255                                r[*itr].insert(i);
256                }
257                return ok;
258        }
259        // ----------------------------------------------------------------------
260        // reverse Hessian sparsity routine called by CppAD
261        bool solve_ode_rev_hes_sparse(
262                size_t                               id ,
263                size_t                                n ,
264                size_t                                m ,
265                size_t                                p ,
266                const vector< std::set<size_t> >&     r ,
267                const vector<bool>&                   s ,
268                vector<bool>&                         t ,
269                const vector< std::set<size_t> >&     u ,
270                vector< std::set<size_t> >&           v )
271        {       // Can just return false if not use RevSparseHes.
272                assert( id == 0 );
273                assert( n == 3 );
274                assert( m == 2 );
275                bool ok = true;
276                std::set<size_t>::const_iterator itr;
277
278                // compute sparsity pattern for T(x) = S(x) * f'(x)
279                vector< std::set<size_t> > S(1);
280                size_t i, j;
281                S[0].clear();
282                for(i = 0; i < m; i++)
283                        if( s[i] )
284                                S[0].insert(i);
285                t = r_ptr_->RevSparseJac(1, s);
286
287                // compute sparsity pattern for A(x)^T = U(x)^T * f'(x)
288                vector< std::set<size_t> > Ut(p), At(p);
289                for(i = 0; i < m; i++)
290                {       for(itr = u[i].begin(); itr != u[i].end(); itr++)
291                                Ut[*itr].insert(i);
292                }
293                At = r_ptr_->RevSparseJac(p, Ut);
294
295                // compute sparsity pattern for H(x)^T = R^T * (S * F)''(x)
296                vector< std::set<size_t> > R(n), Ht(p);
297                for(j = 0; j < n; j++)
298                        R[j] = r[j];
299                r_ptr_->ForSparseJac(p, R);
300                Ht = r_ptr_->RevSparseHes(p, S);
301
302                // compute sparsity pattern for V(x) = A(x) + H(x)^T
303                for(j = 0; j < n; j++)
304                        v[j].clear();
305                for(i = 0; i < p; i++)
306                {       for(itr = At[i].begin(); itr != At[i].end(); itr++)
307                                v[*itr].insert(i);
308                        for(itr = Ht[i].begin(); itr != Ht[i].end(); itr++)
309                                v[*itr].insert(i);
310                }
311
312                // no longer need the forward mode sparsity pattern
313                // (have to reconstruct them every time)
314                r_ptr_->size_forward_set(0);
315
316                return ok;
317        }
318        // ---------------------------------------------------------------------
319        // Declare the AD<double> routine solve_ode(id, ax, ay)
320        CPPAD_USER_ATOMIC(
321                solve_ode                 ,
322                CppAD::vector             ,
323                double                    ,
324                solve_ode_forward         ,
325                solve_ode_reverse         ,
326                solve_ode_for_jac_sparse  ,
327                solve_ode_rev_jac_sparse  ,
328                solve_ode_rev_hes_sparse
329        )
330} // End empty namespace
331
332bool old_usead_2(void)
333{       bool ok = true;
334        using CppAD::NearEqual;
335        double eps = 10. * CppAD::numeric_limits<double>::epsilon();
336
337        // --------------------------------------------------------------------
338        // Create the ADFun<doulbe> r_
339        create_r();
340
341        // --------------------------------------------------------------------
342        // domain and range space vectors
343        size_t n = 3, m = 2;
344        vector< AD<double> > au(n), ax(n), ay(m);
345        au[0]         = 0.0;        // value of z_0 (t) = t, at t = 0
346        ax[1]         = 0.0;        // value of z_1 (t) = t^2/2, at t = 0
347        au[2]         = 1.0;        // final t
348        CppAD::Independent(au);
349        size_t M      = 2;          // number of r steps to take
350        ax[0]         = au[0];      // value of z_0 (t) = t, at t = 0
351        ax[1]         = au[1];      // value of z_1 (t) = t^2/2, at t = 0
352        AD<double> dt = au[2] / M;  // size of each r step
353        ax[2]         = dt;
354        for(size_t i_step = 0; i_step < M; i_step++)
355        {       size_t id = 0;               // not used
356                solve_ode(id, ax, ay);
357                ax[0] = ay[0];
358                ax[1] = ay[1];
359        }
360
361        // create f: u -> y and stop tape recording
362        // y_0(t) = u_0 + t                   = u_0 + u_2
363        // y_1(t) = u_1 + u_0 * t + t^2 / 2   = u_1 + u_0 * u_2 + u_2^2 / 2
364        // where t = u_2
365        ADFun<double> f;
366        f.Dependent(au, ay);
367
368        // --------------------------------------------------------------------
369        // Check forward mode results
370        //
371        // zero order forward
372        vector<double> up(n), yp(m);
373        size_t q  = 0;
374        double u0 = 0.5;
375        double u1 = 0.25;
376        double u2 = 0.75;
377        double check;
378        up[0]     = u0;
379        up[1]     = u1;
380        up[2]     = u2;
381        yp        = f.Forward(q, up);
382        check     = u0 + u2;
383        ok       &= NearEqual( yp[0], check,  eps, eps);
384        check     = u1 + u0 * u2 + u2 * u2 / 2.0;
385        ok       &= NearEqual( yp[1], check,  eps, eps);
386        //
387        // forward mode first derivative w.r.t t
388        q         = 1;
389        up[0]     = 0.0;
390        up[1]     = 0.0;
391        up[2]     = 1.0;
392        yp        = f.Forward(q, up);
393        check     = 1.0;
394        ok       &= NearEqual( yp[0], check,  eps, eps);
395        check     = u0 + u2;
396        ok       &= NearEqual( yp[1], check,  eps, eps);
397        //
398        // forward mode second order Taylor coefficient w.r.t t
399        q         = 2;
400        up[0]     = 0.0;
401        up[1]     = 0.0;
402        up[2]     = 0.0;
403        yp        = f.Forward(q, up);
404        check     = 0.0;
405        ok       &= NearEqual( yp[0], check,  eps, eps);
406        check     = 1.0 / 2.0;
407        ok       &= NearEqual( yp[1], check,  eps, eps);
408        // --------------------------------------------------------------------
409        // reverse mode derivatives of \partial_t y_1 (t)
410        vector<double> w(m * q), dw(n * q);
411        w[0 * q + 0]  = 0.0;
412        w[1 * q + 0]  = 0.0;
413        w[0 * q + 1]  = 0.0;
414        w[1 * q + 1]  = 1.0;
415        dw        = f.Reverse(q, w);
416        // derivative of y_1(u) = u_1 + u_0 * u_2 + u_2^2 / 2,  w.r.t. u
417        // is equal deritative of \partial_u2 y_1(u) w.r.t \partial_u2 u
418        check     = u2;
419        ok       &= NearEqual( dw[0 * q + 1], check,  eps, eps);
420        check     = 1.0;
421        ok       &= NearEqual( dw[1 * q + 1], check,  eps, eps);
422        check     = u0 + u2;
423        ok       &= NearEqual( dw[2 * q + 1], check,  eps, eps);
424        // derivative of \partial_t y_1 w.r.t u = u_0 + t,  w.r.t u
425        check     = 1.0;
426        ok       &= NearEqual( dw[0 * q + 0], check,  eps, eps);
427        check     = 0.0;
428        ok       &= NearEqual( dw[1 * q + 0], check,  eps, eps);
429        check     = 1.0;
430        ok       &= NearEqual( dw[2 * q + 0], check,  eps, eps);
431        // --------------------------------------------------------------------
432        // forward mode sparsity pattern for the Jacobian
433        // f_u = [   1, 0,   1 ]
434        //       [ u_2, 1, u_2 ]
435        size_t i, j, p = n;
436        CppAD::vectorBool r(n * p), s(m * p);
437        // r = identity sparsity pattern
438        for(i = 0; i < n; i++)
439                for(j = 0; j < p; j++)
440                        r[i*n +j] = (i == j);
441        s   = f.ForSparseJac(p, r);
442        ok &= s[ 0 * p + 0] == true;
443        ok &= s[ 0 * p + 1] == false;
444        ok &= s[ 0 * p + 2] == true;
445        ok &= s[ 1 * p + 0] == true;
446        ok &= s[ 1 * p + 1] == true;
447        ok &= s[ 1 * p + 2] == true;
448        // --------------------------------------------------------------------
449        // reverse mode sparsity pattern for the Jacobian
450        q = m;
451        s.resize(q * m);
452        r.resize(q * n);
453        // s = identity sparsity pattern
454        for(i = 0; i < q; i++)
455                for(j = 0; j < m; j++)
456                        s[i*m +j] = (i == j);
457        r   = f.RevSparseJac(q, s);
458        ok &= r[ 0 * n + 0] == true;
459        ok &= r[ 0 * n + 1] == false;
460        ok &= r[ 0 * n + 2] == true;
461        ok &= r[ 1 * n + 0] == true;
462        ok &= r[ 1 * n + 1] == true;
463        ok &= r[ 1 * n + 2] == true;
464
465        // --------------------------------------------------------------------
466        // Hessian sparsity for y_1 (u) = u_1 + u_0 * u_2 + u_2^2 / 2
467        s.resize(m);
468        s[0] = false;
469        s[1] = true;
470        r.resize(n * n);
471        for(i = 0; i < n; i++)
472                for(j = 0; j < n; j++)
473                        r[ i * n + j ] = (i == j);
474        CppAD::vectorBool h(n * n);
475        h   = f.RevSparseHes(n, s);
476        ok &= h[0 * n + 0] == false;
477        ok &= h[0 * n + 1] == false;
478        ok &= h[0 * n + 2] == true;
479        ok &= h[1 * n + 0] == false;
480        ok &= h[1 * n + 1] == false;
481        ok &= h[1 * n + 2] == false;
482        ok &= h[2 * n + 0] == true;
483        ok &= h[2 * n + 1] == false;
484        ok &= h[2 * n + 2] == true;
485
486        // --------------------------------------------------------------------
487        destroy_r();
488
489        // Free all temporary work space associated with old_atomic objects.
490        // (If there are future calls to user atomic functions, they will
491        // create new temporary work space.)
492        CppAD::user_atomic<double>::clear();
493
494        return ok;
495}
496// END C++
Note: See TracBrowser for help on using the repository browser.