source: trunk/example/atomic/reciprocal.cpp @ 2903

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

Add lines from atomic_base documentation to corresponding examples.

package.sh: ignore the bug/build directory.

  • Property svn:keywords set to Id
File size: 13.7 KB
Line 
1// $Id: reciprocal.cpp 2903 2013-09-19 15:53:13Z 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/*
14$begin atomic_reciprocal.cpp$$
15
16$section Reciprocal as an Atomic Operation: Example and Test$$
17$index reciprocal, atomic operation$$
18$index simple, atomic operation$$
19$index atomic, simple operation$$
20$index operation, simple atomic$$
21
22$head Theory$$
23This example demonstrates using $cref atomic_base$$
24to define the operation
25$latex f : \B{R}^n \rightarrow \B{R}^m$$ where
26$latex n = 1$$, $latex m = 1$$, and $latex f(x) = 1 / x$$.
27
28$nospell
29
30$head Start Class Definition$$
31$codep */
32# include <cppad/cppad.hpp>
33namespace {           // isolate items below to this file
34using CppAD::vector;  // abbreviate as vector
35//
36// a utility to compute the union of two sets.
37void my_union(
38        std::set<size_t>&         result  ,
39        const std::set<size_t>&   left    ,
40        const std::set<size_t>&   right   )
41{       std::set<size_t> temp;
42        std::set_union(
43                left.begin()              ,
44                left.end()                ,
45                right.begin()             ,
46                right.end()               ,
47                std::inserter(temp, temp.begin())
48        );
49        result.swap(temp);
50}
51//
52class atomic_reciprocal : public CppAD::atomic_base<double> {
53/* $$
54$head Constructor $$
55$codep */
56        public:
57        // constructor (could use const char* for name)
58        atomic_reciprocal(const std::string& name) : 
59        CppAD::atomic_base<double>(name)
60        { }
61        private:
62/* $$
63$head forward$$
64$codep */
65        // forward mode routine called by CppAD
66        virtual bool forward(
67                size_t                    q ,
68                size_t                    p ,
69                const vector<bool>&      vx ,
70                      vector<bool>&      vy ,
71                const vector<double>&    tx ,
72                      vector<double>&    ty
73        )
74        {       size_t n = tx.size() / (p + 1);
75                size_t m = ty.size() / (p + 1);
76                assert( n == 1 );
77                assert( m == 1 );
78                assert( q <= p );
79
80                // return flag
81                bool ok = p <= 2;
82
83                // check for defining variable information
84                // This case must always be implemented
85                if( vx.size() > 0 )
86                        vy[0] = vx[0];
87
88                // Order zero forward mode.
89                // This case must always be implemented
90                // y^0 = f( x^0 ) = 1 / x^0
91                double f = 1. / tx[0];
92                if( q <= 0 )
93                        ty[0] = f;
94                if( p <= 0 )
95                        return ok;
96                assert( vx.size() == 0 );
97
98                // Order one forward mode.
99                // This case needed if first order forward mode is used.
100                // y^1 = f'( x^0 ) x^1
101                double fp = - f / tx[0];
102                if( q <= 1 )
103                        ty[1] = fp * tx[1]; 
104                if( p <= 1 )
105                        return ok;
106
107                // Order two forward mode.
108                // This case needed if second order forward mode is used.
109                // Y''(t) = X'(t)^\R{T} f''[X(t)] X'(t) + f'[X(t)] X''(t)
110                // 2 y^2  = x^1 * f''( x^0 ) x^1 + 2 f'( x^0 ) x^2
111                double fpp  = - 2.0 * fp / tx[0];
112                ty[2] = tx[1] * fpp * tx[1] / 2.0 + fp * tx[2];
113                if( p <= 2 )
114                        return ok;
115
116                // Assume we are not using forward mode with order > 2
117                assert( ! ok );
118                return ok;
119        }
120/* $$
121$head reverse$$
122$codep */
123        // reverse mode routine called by CppAD
124        virtual bool reverse(
125                size_t                    p ,
126                const vector<double>&    tx ,
127                const vector<double>&    ty ,
128                      vector<double>&    px ,
129                const vector<double>&    py
130        )
131        {       size_t n = tx.size() / (p + 1);
132                size_t m = ty.size() / (p + 1); 
133                assert( px.size() == n * (p + 1) );
134                assert( py.size() == m * (p + 1) );
135                assert( n == 1 );
136                assert( m == 1 );
137                bool ok = p <= 2;       
138
139                double f, fp, fpp, fppp;
140                switch(p)
141                {       case 0:
142                        // This case needed if first order reverse mode is used
143                        // reverse: F^0 ( tx ) = y^0 = f( x^0 )
144                        f     = ty[0];
145                        fp    = - f / tx[0];
146                        px[0] = py[0] * fp;;
147                        assert(ok);
148                        break;
149
150                        case 1:
151                        // This case needed if second order reverse mode is used
152                        // reverse: F^1 ( tx ) = y^1 = f'( x^0 ) x^1
153                        f      = ty[0];
154                        fp     = - f / tx[0];
155                        fpp    = - 2.0 * fp / tx[0];
156                        px[1]  = py[1] * fp;
157                        px[0]  = py[1] * fpp * tx[1];
158                        // reverse: F^0 ( tx ) = y^0 = f( x^0 );
159                        px[0] += py[0] * fp;
160                        assert(ok);
161                        break;
162
163                        case 2:
164                        // This needed if third order reverse mode is used
165                        // reverse: F^2 ( tx ) = y^2 =
166                        //          = x^1 * f''( x^0 ) x^1 / 2 + f'( x^0 ) x^2
167                        f      = ty[0];
168                        fp     = - f / tx[0];
169                        fpp    = - 2.0 * fp / tx[0];
170                        fppp   = - 3.0 * fpp / tx[0];
171                        px[2]  = py[2] * fp;
172                        px[1]  = py[2] * fpp * tx[1];
173                        px[0]  = py[2] * tx[1] * fppp * tx[1] / 2.0 + fpp * tx[2]; 
174                        // reverse: F^1 ( tx ) = y^1 = f'( x^0 ) x^1
175                        px[1] += py[1] * fp;
176                        px[0] += py[1] * fpp * tx[1];
177                        // reverse: F^0 ( tx ) = y^0 = f( x^0 );
178                        px[0] += py[0] * fp;
179                        assert(ok);
180                        break;
181
182                        default:
183                        assert(!ok);
184                }
185                return ok;
186        }
187/* $$
188$head for_sparse_jac$$
189$codep */
190        // forward Jacobian bool sparsity routine called by CppAD
191        virtual bool for_sparse_jac(
192                size_t                                q ,
193                const vector<bool>&                   r ,
194                      vector<bool>&                   s )
195        {       // This function needed if using f.ForSparseJac
196                // with afun.option( CppAD::atomic_base<double>::bool_sparsity_enum )
197                size_t n = r.size() / q;
198                size_t m = s.size() / q;
199                assert( n == 1 );
200                assert( m == 1 );
201
202                // sparsity for S(x) = f'(x) * R is same as sparsity for R
203                for(size_t j = 0; j < q; j++)
204                        s[j] = r[j];
205
206                return true; 
207        }
208        // forward Jacobian set sparsity routine called by CppAD
209        virtual bool for_sparse_jac(
210                size_t                                q ,
211                const vector< std::set<size_t> >&     r ,
212                      vector< std::set<size_t> >&     s )
213        {       // This function needed if using f.ForSparseJac
214                // with afun.option( CppAD::atomic_base<double>::set_sparsity_enum )
215                size_t n = r.size();
216                size_t m = s.size();
217                assert( n == 1 );
218                assert( m == 1 );
219
220                // sparsity for S(x) = f'(x) * R is same as sparsity for R
221                s[0] = r[0];
222
223                return true; 
224        }
225/* $$
226$head rev_sparse_jac$$
227$codep */
228        // reverse Jacobian bool sparsity routine called by CppAD
229        virtual bool rev_sparse_jac(
230                size_t                                q  ,
231                const vector<bool>&                   rt ,
232                      vector<bool>&                   st )
233        {       // This function needed if using RevSparseJac
234                // with afun.option( CppAD::atomic_base<double>::bool_sparsity_enum )
235                size_t n = st.size() / q;
236                size_t m = rt.size() / q;
237                assert( n == 1 );
238                assert( m == 1 );
239
240                // sparsity for S(x)^T = f'(x)^T * R^T is same as sparsity for R^T
241                for(size_t i = 0; i < q; i++)
242                        st[i] = rt[i];
243
244                return true; 
245        }
246        // reverse Jacobian set sparsity routine called by CppAD
247        virtual bool rev_sparse_jac(
248                size_t                                q  ,
249                const vector< std::set<size_t> >&     rt ,
250                      vector< std::set<size_t> >&     st )
251        {       // This function needed if using RevSparseJac
252                // with afun.option( CppAD::atomic_base<double>::set_sparsity_enum )
253                size_t n = st.size();
254                size_t m = rt.size();
255                assert( n == 1 );
256                assert( m == 1 );
257
258                // sparsity for S(x)^T = f'(x)^T * R^T is same as sparsity for R^T
259                st[0] = rt[0];
260
261                return true; 
262        }
263/* $$
264$head rev_sparse_hes$$
265$codep */
266        // reverse Hessian bool sparsity routine called by CppAD
267        virtual bool rev_sparse_hes(
268                const vector<bool>&                   vx,
269                const vector<bool>&                   s ,
270                      vector<bool>&                   t ,
271                size_t                                q ,
272                const vector<bool>&                   r ,
273                const vector<bool>&                   u ,
274                      vector<bool>&                   v )
275        {       // This function needed if using RevSparseHes
276                // with afun.option( CppAD::atomic_base<double>::bool_sparsity_enum )
277                size_t m = s.size();
278                size_t n = t.size();
279                assert( r.size() == n * q );
280                assert( u.size() == m * q );
281                assert( v.size() == n * q );
282                assert( n == 1 );
283                assert( m == 1 );
284
285                // There are no cross term second derivatives for this case,
286                // so it is not necessary to vx.
287
288                // sparsity for T(x) = S(x) * f'(x) is same as sparsity for S
289                t[0] = s[0];
290
291                // V(x) = f'(x)^T * g''(y) * f'(x) * R  +  g'(y) * f''(x) * R
292                // U(x) = g''(y) * f'(x) * R
293                // S(x) = g'(y)
294               
295                // back propagate the sparsity for U, note f'(x) may be non-zero;
296                size_t j;
297                for(j = 0; j < q; j++)
298                        v[j] = u[j];
299
300                // include forward Jacobian sparsity in Hessian sparsity
301                // (note sparsty for f''(x) * R same as for R)
302                if( s[0] )
303                {       for(j = 0; j < q; j++)
304                                v[j] |= r[j];
305                }
306
307                return true;
308        }
309        // reverse Hessian set sparsity routine called by CppAD
310        virtual bool rev_sparse_hes(
311                const vector<bool>&                   vx,
312                const vector<bool>&                   s ,
313                      vector<bool>&                   t ,
314                size_t                                q ,
315                const vector< std::set<size_t> >&     r ,
316                const vector< std::set<size_t> >&     u ,
317                      vector< std::set<size_t> >&     v )
318        {       // This function needed if using RevSparseHes
319                // with afun.option( CppAD::atomic_base<double>::set_sparsity_enum )
320                size_t n = vx.size();
321                size_t m = s.size();
322                assert( t.size() == n );
323                assert( r.size() == n );
324                assert( u.size() == m );
325                assert( v.size() == n );
326                assert( n == 1 );
327                assert( m == 1 );
328
329                // There are no cross term second derivatives for this case,
330                // so it is not necessary to vx.
331
332                // sparsity for T(x) = S(x) * f'(x) is same as sparsity for S
333                t[0] = s[0];
334       
335                // V(x) = f'(x)^T * g''(y) * f'(x) * R  +  g'(y) * f''(x) * R
336                // U(x) = g''(y) * f'(x) * R
337                // S(x) = g'(y)
338               
339                // back propagate the sparsity for U, note f'(x) may be non-zero;
340                v[0] = u[0];
341
342                // include forward Jacobian sparsity in Hessian sparsity
343                // (note sparsty for f''(x) * R same as for R)
344                if( s[0] )
345                        my_union(v[0], v[0], r[0] );
346
347                return true;
348        }
349/* $$
350$head End Class Definition$$
351$codep */
352}; // End of atomic_reciprocal class
353}  // End empty namespace
354
355/* $$
356$head Use Atomic Function$$
357$codep */
358bool reciprocal(void)
359{       bool ok = true;
360        using CppAD::AD;
361        using CppAD::NearEqual;
362        double eps = 10. * CppAD::numeric_limits<double>::epsilon();
363
364        // --------------------------------------------------------------------
365        // Create the atomic reciprocal object
366        atomic_reciprocal afun("atomic_reciprocal");
367        // --------------------------------------------------------------------
368        // Create the function f(x)
369        //
370        // domain space vector
371        size_t n  = 1;
372        double  x0 = 0.5;
373        vector< AD<double> > ax(n);
374        ax[0]     = x0;
375
376        // declare independent variables and start tape recording
377        CppAD::Independent(ax);
378
379        // range space vector
380        size_t m = 1;
381        vector< AD<double> > ay(m);
382
383        // call user function and store reciprocal(x) in au[0]
384        vector< AD<double> > au(m);
385        afun(ax, au);        // u = 1 / x
386
387        // now use AD division to invert to invert the operation
388        ay[0] = 1.0 / au[0]; // y = 1 / u = x
389
390        // create f: x -> y and stop tape recording
391        CppAD::ADFun<double> f;
392        f.Dependent (ax, ay);  // f(x) = x
393
394        // --------------------------------------------------------------------
395        // Check forward mode results
396        //
397        // check function value
398        double check = x0;
399        ok &= NearEqual( Value(ay[0]) , check,  eps, eps);
400
401        // check zero order forward mode
402        size_t p;
403        vector<double> x_p(n), y_p(m);
404        p      = 0;
405        x_p[0] = x0;
406        y_p    = f.Forward(p, x_p);
407        ok &= NearEqual(y_p[0] , check,  eps, eps);
408
409        // check first order forward mode
410        p      = 1;
411        x_p[0] = 1;
412        y_p    = f.Forward(p, x_p);
413        check  = 1.;
414        ok &= NearEqual(y_p[0] , check,  eps, eps);
415
416        // check second order forward mode
417        p      = 2;
418        x_p[0] = 0;
419        y_p    = f.Forward(p, x_p);
420        check  = 0.;
421        ok &= NearEqual(y_p[0] , check,  eps, eps);
422
423        // --------------------------------------------------------------------
424        // Check reverse mode results
425        //
426        // third order reverse mode
427        p     = 3;
428        vector<double> w(m), dw(n * p);
429        w[0]  = 1.;
430        dw    = f.Reverse(p, w);
431        check = 1.;
432        ok &= NearEqual(dw[0] , check,  eps, eps);
433        check = 0.;
434        ok &= NearEqual(dw[1] , check,  eps, eps);
435        ok &= NearEqual(dw[2] , check,  eps, eps);
436
437        // --------------------------------------------------------------------
438        // forward mode sparstiy pattern
439        size_t q = n;
440        CppAD::vectorBool r1(n * q), s1(m * q);
441        r1[0] = true;          // compute sparsity pattern for x[0]
442        //
443        afun.option( CppAD::atomic_base<double>::bool_sparsity_enum );
444        s1    = f.ForSparseJac(q, r1);
445        ok  &= s1[0] == true;  // f[0] depends on x[0] 
446        //
447        afun.option( CppAD::atomic_base<double>::set_sparsity_enum );
448        s1    = f.ForSparseJac(q, r1);
449        ok  &= s1[0] == true;  // f[0] depends on x[0] 
450
451        // --------------------------------------------------------------------
452        // reverse mode sparstiy pattern
453        p = m;
454        CppAD::vectorBool s2(p * m), r2(p * n);
455        s2[0] = true;          // compute sparsity pattern for f[0]
456        //
457        afun.option( CppAD::atomic_base<double>::bool_sparsity_enum );
458        r2    = f.RevSparseJac(p, s2);
459        ok  &= r2[0] == true;  // f[0] depends on x[0] 
460        //
461        afun.option( CppAD::atomic_base<double>::set_sparsity_enum );
462        r2    = f.RevSparseJac(p, s2);
463        ok  &= r2[0] == true;  // f[0] depends on x[0] 
464
465        // --------------------------------------------------------------------
466        // Hessian sparsity (using previous ForSparseJac call)
467        CppAD::vectorBool s3(m), h(q * n);
468        s3[0] = true;        // compute sparsity pattern for f[0]
469        //
470        afun.option( CppAD::atomic_base<double>::bool_sparsity_enum );
471        h     = f.RevSparseHes(q, s3);
472        ok  &= h[0] == true; // second partial of f[0] w.r.t. x[0] may be non-zero
473        //
474        afun.option( CppAD::atomic_base<double>::set_sparsity_enum );
475        h     = f.RevSparseHes(q, s3);
476        ok  &= h[0] == true; // second partial of f[0] w.r.t. x[0] may be non-zero
477
478        // -----------------------------------------------------------------
479        // Free all temporary work space associated with atomic_base objects.
480        // (If there are future calls to user atomic functions, they will
481        // create new temporary work space.)
482        CppAD::atomic_base<double>::clear();
483
484        return ok;
485}
486/* $$
487$$ $comment end nospell$$
488$end
489*/
Note: See TracBrowser for help on using the repository browser.