source: trunk/test_more/old_reciprocal.cpp @ 3788

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

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

commit 431e0a227dbfe6172d265d9d79a2b5b258c5bc18
Author: Brad Bell <bradbell@…>
Date: Tue Feb 9 07:41:29 2016 -0700

  1. Change package.sh to automatically update version on (master branch only).
  2. Change version.sh copy to check and change.

commit a14455414810cfe3c3e4bca90090defc2528a353
Author: Brad Bell <bradbell@…>
Date: Tue Feb 9 06:19:54 2016 -0700

Change check_verbatim to check_srcfile
because all but one of the verbatim commands were changed to srcfile commands.


check_include_omh.sh: old check that file names did not change case (for cygwin development).

commit 4ce45b796b57629332ab46d8ae6df94e0a1ed998
Author: Brad Bell <bradbell@…>
Date: Tue Feb 9 06:04:57 2016 -0700

batch_edit.sh to change $verbatim and $codep to $srcfile and $srccode.


det_by_minor.c: remove some invisible white space.

commit 56553b88c9623c30d2222425a9640b95ce4c8281
Author: Brad Bell <bradbell@…>
Date: Mon Feb 8 18:01:49 2016 -0700

check_jenkins.sh: jenkins.sh no longer takes an argument.
jenkins.sh: fix name of script in error message.

commit 3b8a208cfc7e8ef3c928c17eb291aa3b90ff0050
Author: Brad Bell <bradbell@…>
Date: Mon Feb 8 07:57:02 2016 -0700

new_release.sh: track branches in comments, back to master at OK end, first check of response.

commit 442b7cbc45c022776e8257d3c3404dccdd06c420
Author: Brad Bell <bradbell@…>
Date: Mon Feb 8 06:01:11 2016 -0700

  1. Advance to release 20160000.1.
  2. Check using master version of new_release.sh.
  3. Make sure auto-tools version up to date.
  4. Ask user if doing further testing before commiting new release.

commit f7bdd1f48e72feb05d604da63914022809f45c28
Author: Brad Bell <bradbell@…>
Date: Sun Feb 7 07:59:41 2016 -0700

Add shared library version number to cppad_lib; i.e., cppad_lib.yyyy.mmdd.rel

commit a4c716552e3ad05b337aea58b643c9ad1cbcd4ac
Author: Brad Bell <bradbell@…>
Date: Sun Feb 7 05:25:39 2016 -0700

Make cppad_lib libarary conditional on colpack_prefix being specified.

commit 5e8890eb8de8b0cde146a6ed59c391d7c355ff24
Author: Brad Bell <bradbell@…>
Date: Tue Jan 26 10:49:37 2016 -0700

vector.hpp: fix The -> This.

commit e4e5442b069d7b00e197c31616da32eee20460b3
Merge: c4559d5 ed28b89
Author: Brad Bell <bradbell@…>
Date: Tue Jan 26 09:47:58 2016 -0700

Merge pull request #14 from barak/master


fix spelling in description of cppad_profile_flag on cmake command line

commit ed28b899c9fedab52a578aa7dd73818638081fe6
Author: Barak A. Pearlmutter <barak+git@…>
Date: Tue Jan 26 16:24:32 2016 +0000

typo

File size: 10.6 KB
Line 
1// $Id$
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/*
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$srcfile%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.