# source:trunk/example/atomic/reciprocal.cpp@2899

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

atomic_err_msg.sh: Demonstrate problem with atomic_base error messaging.
get_started.cpp: fix sero->zero.
old_reciprocal.cpp: fix sero->zero.
reciprocal.cpp: fix sero->zero.

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