source: trunk/example/atomic/tangent.cpp @ 2861

Last change on this file since 2861 was 2861, checked in by bradbell, 7 years ago
  1. Use namespace in call to base constructor (works on older compilers).
  2. Make destruction virtual (avoid warning on some compilers).
  • Property svn:keywords set to Id
File size: 13.4 KB
Line 
1// $Id: tangent.cpp 2861 2013-05-28 20:29:31Z 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$code
30$verbatim%example/atomic/tangent.cpp%0%// BEGIN C++%// END C++%1%$$
31$$
32
33$end
34*/
35// BEGIN C++
36# include <cppad/cppad.hpp>
37
38namespace { // Begin empty namespace
39using CppAD::vector;
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_tangent : public CppAD::atomic_base<float> {
58private:
59        const bool hyperbolic_; // is this hyperbolic tangent
60public:
61        // ----------------------------------------------------------------------
62        // constructor
63        atomic_tangent(const char* name, bool hyperbolic) 
64        : CppAD::atomic_base<float>(name),
65        hyperbolic_(hyperbolic)
66        { }
67private:
68        // ----------------------------------------------------------------------
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        // reverse mode routine called by CppAD
125        virtual bool reverse(
126                size_t                    p ,
127                const vector<float>&     tx ,
128                const vector<float>&    tzy ,
129                      vector<float>&     px ,
130                const vector<float>&    pzy
131        )
132        {       size_t p1 = p + 1;
133                size_t n  = tx.size()  / p1;
134                size_t m  = tzy.size() / p1;   
135                assert( px.size()  == n * p1 );
136                assert( pzy.size() == m * p1 );
137                assert( n == 1 );
138                assert( m == 2 );
139
140                size_t j, k;
141
142                // copy because partials w.r.t. y and z need to change
143                vector<float> qzy = pzy;
144
145                // initialize accumultion of reverse mode partials
146                for(k = 0; k < p1; k++)
147                        px[k] = 0.;
148
149                // eliminate positive orders
150                for(j = p; j > 0; j--)
151                {       float j_inv = 1.f / float(j);
152                        if( hyperbolic_ )
153                                j_inv = - j_inv;
154
155                        // H_{x^{(k)}} += delta(j-k) +- H_{z^{(j)} y^{(j-k)} * k / j
156                        px[j] += qzy[j];
157                        for(k = 1; k <= j; k++)
158                                px[k] += qzy[j] * tzy[p1 + j-k] * k * j_inv; 
159
160                        // H_{y^{j-k)} += +- H_{z^{(j)} x^{(k)} * k / j
161                        for(k = 1; k <= j; k++)
162                                qzy[p1 + j-k] += qzy[j] * tx[k] * k * j_inv; 
163
164                        // H_{z^{(k)}} += H_{y^{(j-1)}} * z^{(j-k-1)} * 2.
165                        for(k = 0; k < j; k++)
166                                qzy[k] += qzy[p1 + j-1] * tzy[j-k-1] * 2.f; 
167                }
168
169                // eliminate order zero
170                if( hyperbolic_ )
171                        px[0] += qzy[0] * (1.f - tzy[p1 + 0]);
172                else
173                        px[0] += qzy[0] * (1.f + tzy[p1 + 0]);
174
175                return true; 
176        }
177        // ----------------------------------------------------------------------
178        // forward Jacobian sparsity routine called by CppAD
179        virtual bool for_sparse_jac(
180                size_t                                q ,
181                const vector<bool>&                   r ,
182                      vector<bool>&                   s )
183        {       size_t n = r.size() / q;
184                size_t m = s.size() / q;
185                assert( n == 1 );
186                assert( m == 2 );
187
188                // sparsity for S(x) = f'(x) * R
189                for(size_t j = 0; j < q; j++)
190                {       s[0 * q + j] = r[j];
191                        s[1 * q + j] = r[j];
192                }
193
194                return true;
195        }
196        // forward Jacobian sparsity routine called by CppAD
197        virtual bool for_sparse_jac(
198                size_t                                q ,
199                const vector< std::set<size_t> >&     r ,
200                      vector< std::set<size_t> >&     s )
201        {       size_t n = r.size();
202                size_t m = s.size();
203                assert( n == 1 );
204                assert( m == 2 );
205
206                // sparsity for S(x) = f'(x) * R
207                s[0] = r[0];
208                s[1] = r[0];
209
210                return true;
211        }
212        // ----------------------------------------------------------------------
213        // reverse Jacobian sparsity routine called by CppAD
214        virtual bool rev_sparse_jac(
215                size_t                                q ,
216                const vector<bool>&                  rt ,
217                      vector<bool>&                  st )
218        {       size_t n = st.size() / q;
219                size_t m = rt.size() / q;
220                assert( n == 1 );
221                assert( m == 2 );
222
223                // sparsity for S(x)^T = f'(x)^T * R^T
224                for(size_t j = 0; j < q; j++)
225                        st[j] = rt[0 * q + j] | rt[1 * q + j];
226
227                return true; 
228        }
229        // reverse Jacobian sparsity routine called by CppAD
230        virtual bool rev_sparse_jac(
231                size_t                                q ,
232                const vector< std::set<size_t> >&    rt ,
233                      vector< std::set<size_t> >&    st )
234        {       size_t n = st.size();
235                size_t m = rt.size();
236                assert( n == 1 );
237                assert( m == 2 );
238
239                // sparsity for S(x)^T = f'(x)^T * R^T
240                my_union(st[0], rt[0], rt[1]);
241                return true; 
242        }
243        // ----------------------------------------------------------------------
244        // reverse Hessian sparsity routine called by CppAD
245        virtual bool rev_sparse_hes(
246                const vector<bool>&                   vx,
247                const vector<bool>&                   s ,
248                      vector<bool>&                   t ,
249                size_t                                q ,
250                const vector<bool>&                   r ,
251                const vector<bool>&                   u ,
252                      vector<bool>&                   v )
253        {
254                size_t m = s.size();
255                size_t n = t.size();
256                assert( r.size() == n * q );
257                assert( u.size() == m * q );
258                assert( v.size() == n * q );
259                assert( n == 1 );
260                assert( m == 2 );
261
262                // There are no cross term second derivatives for this case,
263                // so it is not necessary to vx.
264
265                // sparsity for T(x) = S(x) * f'(x)
266                t[0] =  s[0] | s[1];
267
268                // V(x) = f'(x)^T * g''(y) * f'(x) * R  +  g'(y) * f''(x) * R
269                // U(x) = g''(y) * f'(x) * R
270                // S(x) = g'(y)
271               
272                // back propagate the sparsity for U, note both components
273                // of f'(x) may be non-zero;
274                size_t j;
275                for(j = 0; j < q; j++)
276                        v[j] = u[ 0 * q + j ] | u[ 1 * q + j ];
277
278                // include forward Jacobian sparsity in Hessian sparsity
279                // (note sparsty for f''(x) * R same as for R)
280                if( s[0] | s[1] )
281                {       for(j = 0; j < q; j++)
282                                v[j] |= r[j];
283                }
284
285                return true;
286        }
287        // reverse Hessian sparsity routine called by CppAD
288        virtual bool rev_sparse_hes(
289                const vector<bool>&                   vx,
290                const vector<bool>&                   s ,
291                      vector<bool>&                   t ,
292                size_t                                q ,
293                const vector< std::set<size_t> >&     r ,
294                const vector< std::set<size_t> >&     u ,
295                      vector< std::set<size_t> >&     v )
296        {       size_t m = s.size();
297                size_t n = t.size();
298                assert( r.size() == n );
299                assert( u.size() == m );
300                assert( v.size() == n );
301                assert( n == 1 );
302                assert( m == 2 );
303
304                // There are no cross term second derivatives for this case,
305                // so it is not necessary to vx.
306
307                // sparsity for T(x) = S(x) * f'(x)
308                t[0] =  s[0] | s[1];
309
310                // V(x) = f'(x)^T * g''(y) * f'(x) * R  +  g'(y) * f''(x) * R
311                // U(x) = g''(y) * f'(x) * R
312                // S(x) = g'(y)
313               
314                // back propagate the sparsity for U, note both components
315                // of f'(x) may be non-zero;
316                my_union(v[0], u[0], u[1]);
317
318                // include forward Jacobian sparsity in Hessian sparsity
319                // (note sparsty for f''(x) * R same as for R)
320                if( s[0] | s[1] )
321                        my_union(v[0], v[0], r[0]);
322
323                return true;
324        }
325}; // End of atomic_tangent class
326}  // End empty namespace
327
328bool tangent(void)
329{       bool ok = true;
330        using CppAD::AD;
331        using CppAD::NearEqual;
332        float eps = 10.f * CppAD::numeric_limits<float>::epsilon();
333
334        // --------------------------------------------------------------------
335        // Creater a tan and tanh object
336        atomic_tangent my_tan("my_tan", false), my_tanh("my_tanh", true);
337        // --------------------------------------------------------------------
338
339        // domain space vector
340        size_t n  = 1;
341        float  x0 = 0.5;
342        CppAD::vector< AD<float> > ax(n);
343        ax[0]     = x0;
344
345        // declare independent variables and start tape recording
346        CppAD::Independent(ax);
347
348        // range space vector
349        size_t m = 3;
350        CppAD::vector< AD<float> > af(m);
351
352        // temporary vector for computations
353        // (my_tan and my_tanh computes tan or tanh and its square)
354        CppAD::vector< AD<float> > az(2);
355
356        // call atomic tan function and store tan(x) in f[0] (ignore tan(x)^2)
357        my_tan(ax, az);
358        af[0] = az[0];
359
360        // call atomic tanh function and store tanh(x) in f[1] (ignore tanh(x)^2)
361        my_tanh(ax, az);
362        af[1] = az[0];
363
364        // put a constant in f[2] = tanh(1.) (for sparsity pattern testing)
365        CppAD::vector< AD<float> > one(1);
366        one[0] = 1.;
367        my_tanh(one, az);
368        af[2] = az[0]; 
369
370        // create f: x -> f and stop tape recording
371        CppAD::ADFun<float> F;
372        F.Dependent(ax, af); 
373
374        // check function value
375        float tan = std::tan(x0);
376        ok &= NearEqual(af[0] , tan,  eps, eps);
377        float tanh = std::tanh(x0);
378        ok &= NearEqual(af[1] , tanh,  eps, eps);
379
380        // check zero order forward
381        CppAD::vector<float> x(n), f(m);
382        x[0] = x0;
383        f    = F.Forward(0, x);
384        ok &= NearEqual(f[0] , tan,  eps, eps);
385        ok &= NearEqual(f[1] , tanh,  eps, eps);
386
387        // compute first partial of f w.r.t. x[0] using forward mode
388        CppAD::vector<float> dx(n), df(m);
389        dx[0] = 1.;
390        df    = F.Forward(1, dx);
391
392        // compute derivative of tan - tanh using reverse mode
393        CppAD::vector<float> w(m), dw(n);
394        w[0]  = 1.;
395        w[1]  = 1.;
396        w[2]  = 0.;
397        dw    = F.Reverse(1, w);
398
399        // tan'(x)   = 1 + tan(x)  * tan(x)
400        // tanh'(x)  = 1 - tanh(x) * tanh(x)
401        float tanp  = 1.f + tan * tan; 
402        float tanhp = 1.f - tanh * tanh; 
403        ok   &= NearEqual(df[0], tanp, eps, eps);
404        ok   &= NearEqual(df[1], tanhp, eps, eps);
405        ok   &= NearEqual(dw[0], w[0]*tanp + w[1]*tanhp, eps, eps);
406
407        // compute second partial of f w.r.t. x[0] using forward mode
408        CppAD::vector<float> ddx(n), ddf(m);
409        ddx[0] = 0.;
410        ddf    = F.Forward(2, ddx);
411
412        // compute second derivative of tan - tanh using reverse mode
413        CppAD::vector<float> ddw(2);
414        ddw   = F.Reverse(2, w);
415
416        // tan''(x)   = 2 *  tan(x) * tan'(x)
417        // tanh''(x)  = - 2 * tanh(x) * tanh'(x)
418        // Note that second order Taylor coefficient for u half the
419        // corresponding second derivative.
420        float two    = 2;
421        float tanpp  =   two * tan * tanp;
422        float tanhpp = - two * tanh * tanhp;
423        ok   &= NearEqual(two * ddf[0], tanpp, eps, eps);
424        ok   &= NearEqual(two * ddf[1], tanhpp, eps, eps);
425        ok   &= NearEqual(ddw[0], w[0]*tanp  + w[1]*tanhp , eps, eps);
426        ok   &= NearEqual(ddw[1], w[0]*tanpp + w[1]*tanhpp, eps, eps);
427
428        // Forward mode computation of sparsity pattern for F.
429        size_t q = n;
430        // user vectorBool because m and n are small
431        CppAD::vectorBool r1(q), s1(m * q);
432        r1[0] = true;            // propagate sparsity for x[0]
433        s1    = F.ForSparseJac(q, r1);
434        ok  &= (s1[0] == true);  // f[0] depends on x[0]
435        ok  &= (s1[1] == true);  // f[1] depends on x[0]
436        ok  &= (s1[2] == false); // f[2] does not depend on x[0]
437
438        // Reverse mode computation of sparsity pattern for F.
439        size_t p = m;
440        CppAD::vectorBool s2(p * m), r2(p * n);
441        // Sparsity pattern for identity matrix
442        size_t i, j;
443        for(i = 0; i < p; i++)
444        {       for(j = 0; j < m; j++)
445                        s2[i * p + j] = (i == j);
446        }
447        r2   = F.RevSparseJac(p, s2);
448        ok  &= (r2[0] == true);  // f[0] depends on x[0]
449        ok  &= (r2[1] == true);  // f[1] depends on x[0]
450        ok  &= (r2[2] == false); // f[2] does not depend on x[0]
451
452        // Hessian sparsity for f[0]
453        CppAD::vectorBool s3(m), h(q * n);
454        s3[0] = true;
455        s3[1] = false;
456        s3[2] = false;
457        h    = F.RevSparseHes(q, s3);
458        ok  &= (h[0] == true);  // Hessian is non-zero
459
460        // Hessian sparsity for f[2]
461        s3[0] = false;
462        s3[2] = true;
463        h    = F.RevSparseHes(q, s3);
464        ok  &= (h[0] == false);  // Hessian is zero
465
466        // check tanh results for a large value of x
467        x[0]  = std::numeric_limits<float>::max() / two;
468        f     = F.Forward(0, x);
469        tanh  = 1.;
470        ok   &= NearEqual(f[1], tanh, eps, eps);
471        df    = F.Forward(1, dx);
472        tanhp = 0.;
473        ok   &= NearEqual(df[1], tanhp, eps, eps);
474 
475        // --------------------------------------------------------------------
476        // Free all temporary work space associated with atomic_basen objects.
477        // (If there are future calls to user atomic functions, they will
478        // create new temporary work space.)
479        CppAD::user_atomic<float>::clear();
480
481        return ok;
482}
483// END C++
Note: See TracBrowser for help on using the repository browser.