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