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