source: trunk/test_more/old_tan.cpp @ 3741

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

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

commit cd48f5689e158e2758a3bfb57f2e6285b62b8ec1
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 18:54:45 2015 -0700

Run reduce_index on files that were changed today.


reduce_index.py: add some more words to ignore.

commit cf4a8e3a59f06d919ef9a2d911d7c1e462896b39
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 18:35:02 2015 -0700

Change depreciated to deprecated (typo).

commit 37e09901ab78be253fa50004cdc041b042015eeb
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 18:17:29 2015 -0700

Remove invisible white space.

commit 5ffd5cb372caac6c371af5718b8bfc437c8d7f97
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 18:17:09 2015 -0700

Add the deprecation date to some features that were missing it.

commit de437774c18eda2c619cfd8676985113207eb736
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 15:35:52 2015 -0700

whats_new_15.omh: update from previous commit.
test_more.cpp: alphabetize and white space.

commit e8f32f25fe12edf8a8ffa29eb8d2848fef007bd1
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 15:21:17 2015 -0700

  1. Alphabatize list.
  2. remove old_mat_mul.hpp (deprecated).
  3. Correct $cref -> $rref.

commit f205c22e6c4bcabd516214d40e18c4bc4d61a9d1
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 14:44:00 2015 -0700

Remove invisible white space.

commit 1a73e41a8880f8fe351ba4d620cd25055b2a8e8d
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 14:43:43 2015 -0700

Move some deprecated examples from example directory to test_more directory.

commit 2ea5990965e3f8fd00a6393a60b0a59152b41e15
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 14:36:34 2015 -0700

check_user_def.sh: must not include check_user_def.sh in files searched.

commit 26ae3bc6ec918032701b7c1aa37c5b8f9acc6d55
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 13:16:26 2015 -0700

Remove invisible white space.

commit 6cbc289840bc232aff1f676d9ab3b795028aaa36
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 13:16:22 2015 -0700

Add some of the missing dates when certain features were deprecated.

commit 8e5544f467fec5d60eeb1a3739fc024bd5305014
Author: Brad Bell <bradbell@…>
Date: Tue Oct 6 09:23:13 2015 -0700

Automatic check of user API preprocessor symbols and its indicated corrections.

File size: 11.8 KB
Line 
1// $Id$
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-15 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_tan.cpp$$
15$spell
16        Tanh
17$$
18
19$section Old Tan and Tanh as User Atomic Operations: Example and Test$$
20
21$head Deprecated 2013-05-27$$
22This example has not deprecated;
23see $cref atomic_tangent.cpp$$ instead.
24
25$head Theory$$
26The code below uses the $cref tan_forward$$ and $cref tan_reverse$$
27to implement the tangent ($icode%id% == 0%$$) and hyperbolic tangent
28($icode%id% == 1%$$) functions as user atomic operations.
29
30$code
31$verbatim%test_more/old_tan.cpp%0%// BEGIN C++%// END C++%1%$$
32$$
33
34$end
35*/
36// BEGIN C++
37# include <cppad/cppad.hpp>
38
39namespace { // Begin empty namespace
40        using CppAD::vector;
41
42        // a utility to compute the union of two sets.
43        void my_union(
44                std::set<size_t>&         result  ,
45                const std::set<size_t>&   left    ,
46                const std::set<size_t>&   right   )
47        {       std::set<size_t> temp;
48                std::set_union(
49                        left.begin()              ,
50                        left.end()                ,
51                        right.begin()             ,
52                        right.end()               ,
53                        std::inserter(temp, temp.begin())
54                );
55                result.swap(temp);
56        }
57
58        // ----------------------------------------------------------------------
59        // forward mode routine called by CppAD
60        bool old_tan_forward(
61                size_t                   id ,
62                size_t                order ,
63                size_t                    n ,
64                size_t                    m ,
65                const vector<bool>&      vx ,
66                vector<bool>&           vzy ,
67                const vector<float>&     tx ,
68                vector<float>&          tzy
69        )
70        {
71                assert( id == 0 || id == 1 );
72                assert( n == 1 );
73                assert( m == 2 );
74                assert( tx.size() >= (order+1) * n );
75                assert( tzy.size() >= (order+1) * m );
76
77                size_t n_order = order + 1;
78                size_t j = order;
79                size_t k;
80
81                // check if this is during the call to old_tan(id, ax, ay)
82                if( vx.size() > 0 )
83                {       assert( vx.size() >= n );
84                        assert( vzy.size() >= m );
85
86                        // now setvzy
87                        vzy[0] = vx[0];
88                        vzy[1] = vx[0];
89                }
90
91                if( j == 0 )
92                {       // z^{(0)} = tan( x^{(0)} ) or tanh( x^{(0)} )
93                        if( id == 0 )
94                                tzy[0] = tan( tx[0] );
95                        else    tzy[0] = tanh( tx[0] );
96
97                        // y^{(0)} = z^{(0)} * z^{(0)}
98                        tzy[n_order + 0] = tzy[0] * tzy[0];
99                }
100                else
101                {       float j_inv = 1.f / float(j);
102                        if( id == 1 )
103                                j_inv = - j_inv;
104
105                        // z^{(j)} = x^{(j)} +- sum_{k=1}^j k x^{(k)} y^{(j-k)} / j
106                        tzy[j] = tx[j];
107                        for(k = 1; k <= j; k++)
108                                tzy[j] += tx[k] * tzy[n_order + j-k] * k * j_inv;
109
110                        // y^{(j)} = sum_{k=0}^j z^{(k)} z^{(j-k)}
111                        tzy[n_order + j] = 0.;
112                        for(k = 0; k <= j; k++)
113                                tzy[n_order + j] += tzy[k] * tzy[j-k];
114                }
115
116                // All orders are implemented and there are no possible errors
117                return true;
118        }
119        // ----------------------------------------------------------------------
120        // reverse mode routine called by CppAD
121        bool old_tan_reverse(
122                size_t                   id ,
123                size_t                order ,
124                size_t                    n ,
125                size_t                    m ,
126                const vector<float>&     tx ,
127                const vector<float>&    tzy ,
128                vector<float>&           px ,
129                const vector<float>&    pzy
130        )
131        {
132                assert( id == 0 || id == 1 );
133                assert( n == 1 );
134                assert( m == 2 );
135                assert( tx.size() >= (order+1) * n );
136                assert( tzy.size() >= (order+1) * m );
137                assert( px.size() >= (order+1) * n );
138                assert( pzy.size() >= (order+1) * m );
139
140                size_t n_order = order + 1;
141                size_t j, k;
142
143                // copy because partials w.r.t. y and z need to change
144                vector<float> qzy = pzy;
145
146                // initialize accumultion of reverse mode partials
147                for(k = 0; k < n_order; k++)
148                        px[k] = 0.;
149
150                // eliminate positive orders
151                for(j = order; j > 0; j--)
152                {       float j_inv = 1.f / float(j);
153                        if( id == 1 )
154                                j_inv = - j_inv;
155
156                        // H_{x^{(k)}} += delta(j-k) +- H_{z^{(j)} y^{(j-k)} * k / j
157                        px[j] += qzy[j];
158                        for(k = 1; k <= j; k++)
159                                px[k] += qzy[j] * tzy[n_order + j-k] * k * j_inv;
160
161                        // H_{y^{j-k)} += +- H_{z^{(j)} x^{(k)} * k / j
162                        for(k = 1; k <= j; k++)
163                                qzy[n_order + j-k] += qzy[j] * tx[k] * k * j_inv;
164
165                        // H_{z^{(k)}} += H_{y^{(j-1)}} * z^{(j-k-1)} * 2.
166                        for(k = 0; k < j; k++)
167                                qzy[k] += qzy[n_order + j-1] * tzy[j-k-1] * 2.f;
168                }
169
170                // eliminate order zero
171                if( id == 0 )
172                        px[0] += qzy[0] * (1.f + tzy[n_order + 0]);
173                else
174                        px[0] += qzy[0] * (1.f - tzy[n_order + 0]);
175
176                return true;
177        }
178        // ----------------------------------------------------------------------
179        // forward Jacobian sparsity routine called by CppAD
180        bool old_tan_for_jac_sparse(
181                size_t                               id ,
182                size_t                                n ,
183                size_t                                m ,
184                size_t                                p ,
185                const vector< std::set<size_t> >&     r ,
186                vector< std::set<size_t> >&           s )
187        {
188                assert( n == 1 );
189                assert( m == 2 );
190                assert( id == 0 || id == 1 );
191                assert( r.size() >= n );
192                assert( s.size() >= m );
193
194                // sparsity for z and y are the same as for x
195                s[0] = r[0];
196                s[1] = r[0];
197
198                return true;
199        }
200        // ----------------------------------------------------------------------
201        // reverse Jacobian sparsity routine called by CppAD
202        bool old_tan_rev_jac_sparse(
203                size_t                               id ,
204                size_t                                n ,
205                size_t                                m ,
206                size_t                                p ,
207                vector< std::set<size_t> >&           r ,
208                const vector< std::set<size_t> >&     s )
209        {
210                assert( n == 1 );
211                assert( m == 2 );
212                assert( id == 0 || id == 1 );
213                assert( r.size() >= n );
214                assert( s.size() >= m );
215
216                // note that, if the users code only uses z, and not y,
217                // we could just set r[0] = s[0]
218                my_union(r[0], s[0], s[1]);
219                return true;
220        }
221        // ----------------------------------------------------------------------
222        // reverse Hessian sparsity routine called by CppAD
223        bool old_tan_rev_hes_sparse(
224                size_t                               id ,
225                size_t                                n ,
226                size_t                                m ,
227                size_t                                p ,
228                const vector< std::set<size_t> >&     r ,
229                const vector<bool>&                   s ,
230                vector<bool>&                         t ,
231                const vector< std::set<size_t> >&     u ,
232                vector< std::set<size_t> >&           v )
233        {
234                assert( n == 1 );
235                assert( m == 2 );
236                assert( id == 0 || id == 1 );
237                assert( r.size() >= n );
238                assert( s.size() >= m );
239                assert( t.size() >= n );
240                assert( u.size() >= m );
241                assert( v.size() >= n );
242
243                // back propagate Jacobian sparsity. If users code only uses z,
244                // we could just set t[0] = s[0];
245                t[0] =  s[0] | s[1];
246
247                // back propagate Hessian sparsity, ...
248                my_union(v[0], u[0], u[1]);
249
250                // convert forward Jacobian sparsity to Hessian sparsity
251                // because tan and tanh are nonlinear
252                if( t[0] )
253                        my_union(v[0], v[0], r[0]);
254
255                return true;
256        }
257        // ---------------------------------------------------------------------
258        // Declare the AD<float> routine old_tan(id, ax, ay)
259        CPPAD_USER_ATOMIC(
260                old_tan                 ,
261                CppAD::vector           ,
262                float                   ,
263                old_tan_forward         ,
264                old_tan_reverse         ,
265                old_tan_for_jac_sparse  ,
266                old_tan_rev_jac_sparse  ,
267                old_tan_rev_hes_sparse
268        )
269} // End empty namespace
270
271bool old_tan(void)
272{       bool ok = true;
273        using CppAD::AD;
274        using CppAD::NearEqual;
275        float eps = 10.f * CppAD::numeric_limits<float>::epsilon();
276
277        // domain space vector
278        size_t n  = 1;
279        float  x0 = 0.5;
280        CppAD::vector< AD<float> > ax(n);
281        ax[0]     = x0;
282
283        // declare independent variables and start tape recording
284        CppAD::Independent(ax);
285
286        // range space vector
287        size_t m = 3;
288        CppAD::vector< AD<float> > af(m);
289
290        // temporary vector for old_tan computations
291        // (old_tan computes tan or tanh and its square)
292        CppAD::vector< AD<float> > az(2);
293
294        // call user tan function and store tan(x) in f[0] (ignore tan(x)^2)
295        size_t id = 0;
296        old_tan(id, ax, az);
297        af[0] = az[0];
298
299        // call user tanh function and store tanh(x) in f[1] (ignore tanh(x)^2)
300        id = 1;
301        old_tan(id, ax, az);
302        af[1] = az[0];
303
304        // put a constant in f[2] = tanh(1.) (for sparsity pattern testing)
305        CppAD::vector< AD<float> > one(1);
306        one[0] = 1.;
307        old_tan(id, one, az);
308        af[2] = az[0];
309
310        // create f: x -> f and stop tape recording
311        CppAD::ADFun<float> F;
312        F.Dependent(ax, af);
313
314        // check function value
315        float tan = std::tan(x0);
316        ok &= NearEqual(af[0] , tan,  eps, eps);
317        float tanh = std::tanh(x0);
318        ok &= NearEqual(af[1] , tanh,  eps, eps);
319
320        // check zero order forward
321        CppAD::vector<float> x(n), f(m);
322        x[0] = x0;
323        f    = F.Forward(0, x);
324        ok &= NearEqual(f[0] , tan,  eps, eps);
325        ok &= NearEqual(f[1] , tanh,  eps, eps);
326
327        // compute first partial of f w.r.t. x[0] using forward mode
328        CppAD::vector<float> dx(n), df(m);
329        dx[0] = 1.;
330        df    = F.Forward(1, dx);
331
332        // compute derivative of tan - tanh using reverse mode
333        CppAD::vector<float> w(m), dw(n);
334        w[0]  = 1.;
335        w[1]  = 1.;
336        w[2]  = 0.;
337        dw    = F.Reverse(1, w);
338
339        // tan'(x)   = 1 + tan(x)  * tan(x)
340        // tanh'(x)  = 1 - tanh(x) * tanh(x)
341        float tanp  = 1.f + tan * tan;
342        float tanhp = 1.f - tanh * tanh;
343        ok   &= NearEqual(df[0], tanp, eps, eps);
344        ok   &= NearEqual(df[1], tanhp, eps, eps);
345        ok   &= NearEqual(dw[0], w[0]*tanp + w[1]*tanhp, eps, eps);
346
347        // compute second partial of f w.r.t. x[0] using forward mode
348        CppAD::vector<float> ddx(n), ddf(m);
349        ddx[0] = 0.;
350        ddf    = F.Forward(2, ddx);
351
352        // compute second derivative of tan - tanh using reverse mode
353        CppAD::vector<float> ddw(2);
354        ddw   = F.Reverse(2, w);
355
356        // tan''(x)   = 2 *  tan(x) * tan'(x)
357        // tanh''(x)  = - 2 * tanh(x) * tanh'(x)
358        // Note that second order Taylor coefficient for u half the
359        // corresponding second derivative.
360        float two    = 2;
361        float tanpp  =   two * tan * tanp;
362        float tanhpp = - two * tanh * tanhp;
363        ok   &= NearEqual(two * ddf[0], tanpp, eps, eps);
364        ok   &= NearEqual(two * ddf[1], tanhpp, eps, eps);
365        ok   &= NearEqual(ddw[0], w[0]*tanp  + w[1]*tanhp , eps, eps);
366        ok   &= NearEqual(ddw[1], w[0]*tanpp + w[1]*tanhpp, eps, eps);
367
368        // Forward mode computation of sparsity pattern for F.
369        size_t p = n;
370        // user vectorBool because m and n are small
371        CppAD::vectorBool r1(p), s1(m * p);
372        r1[0] = true;            // propagate sparsity for x[0]
373        s1    = F.ForSparseJac(p, r1);
374        ok  &= (s1[0] == true);  // f[0] depends on x[0]
375        ok  &= (s1[1] == true);  // f[1] depends on x[0]
376        ok  &= (s1[2] == false); // f[2] does not depend on x[0]
377
378        // Reverse mode computation of sparsity pattern for F.
379        size_t q = m;
380        CppAD::vectorBool s2(q * m), r2(q * n);
381        // Sparsity pattern for identity matrix
382        size_t i, j;
383        for(i = 0; i < q; i++)
384        {       for(j = 0; j < m; j++)
385                        s2[i * q + j] = (i == j);
386        }
387        r2   = F.RevSparseJac(q, s2);
388        ok  &= (r2[0] == true);  // f[0] depends on x[0]
389        ok  &= (r2[1] == true);  // f[1] depends on x[0]
390        ok  &= (r2[2] == false); // f[2] does not depend on x[0]
391
392        // Hessian sparsity for f[0]
393        CppAD::vectorBool s3(m), h(p * n);
394        s3[0] = true;
395        s3[1] = false;
396        s3[2] = false;
397        h    = F.RevSparseHes(p, s3);
398        ok  &= (h[0] == true);  // Hessian is non-zero
399
400        // Hessian sparsity for f[2]
401        s3[0] = false;
402        s3[2] = true;
403        h    = F.RevSparseHes(p, s3);
404        ok  &= (h[0] == false);  // Hessian is zero
405
406        // check tanh results for a large value of x
407        x[0]  = std::numeric_limits<float>::max() / two;
408        f     = F.Forward(0, x);
409        tanh  = 1.;
410        ok   &= NearEqual(f[1], tanh, eps, eps);
411        df    = F.Forward(1, dx);
412        tanhp = 0.;
413        ok   &= NearEqual(df[1], tanhp, eps, eps);
414
415        // --------------------------------------------------------------------
416        // Free all temporary work space associated with old_atomic objects.
417        // (If there are future calls to user atomic functions, they will
418        // create new temporary work space.)
419        CppAD::user_atomic<float>::clear();
420
421        return ok;
422}
423// END C++
Note: See TracBrowser for help on using the repository browser.