source: trunk/cppad/local/exp_op.hpp @ 3301

Last change on this file since 3301 was 3301, checked in by bradbell, 6 years ago

merge in multiple forward direcitons from branches/forward_dir

  • Property svn:keywords set to Id
File size: 4.6 KB
Line 
1/* $Id: exp_op.hpp 3301 2014-05-24 05:20:21Z bradbell $ */
2# ifndef CPPAD_EXP_OP_INCLUDED
3# define CPPAD_EXP_OP_INCLUDED
4
5/* --------------------------------------------------------------------------
6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 Bradley M. Bell
7
8CppAD is distributed under multiple licenses. This distribution is under
9the terms of the
10                    Eclipse Public License Version 1.0.
11
12A copy of this license is included in the COPYING file of this distribution.
13Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
14-------------------------------------------------------------------------- */
15
16
17namespace CppAD { // BEGIN_CPPAD_NAMESPACE
18/*!
19\file exp_op.hpp
20Forward and reverse mode calculations for z = exp(x).
21*/
22
23
24/*!
25Forward mode Taylor coefficient for result of op = ExpOp.
26
27The C++ source code corresponding to this operation is
28\verbatim
29        z = exp(x)
30\endverbatim
31
32\copydetails forward_unary1_op
33*/
34template <class Base>
35inline void forward_exp_op(
36        size_t p           ,
37        size_t q           ,
38        size_t i_z         ,
39        size_t i_x         ,
40        size_t cap_order   , 
41        Base*  taylor      )
42{       
43        // check assumptions
44        CPPAD_ASSERT_UNKNOWN( NumArg(ExpOp) == 1 );
45        CPPAD_ASSERT_UNKNOWN( NumRes(ExpOp) == 1 );
46        CPPAD_ASSERT_UNKNOWN( i_x < i_z );
47        CPPAD_ASSERT_UNKNOWN( q < cap_order );
48        CPPAD_ASSERT_UNKNOWN( p <= q );
49
50        // Taylor coefficients corresponding to argument and result
51        Base* x = taylor + i_x * cap_order;
52        Base* z = taylor + i_z * cap_order;
53
54        size_t k;
55        if( p == 0 )
56        {       z[0] = exp( x[0] );
57                p++;
58        }
59        for(size_t j = p; j <= q; j++)
60        {
61                z[j] = x[1] * z[j-1];
62                for(k = 2; k <= j; k++)
63                        z[j] += Base(k) * x[k] * z[j-k];
64                z[j] /= Base(j);
65        }
66}
67
68
69/*!
70Multiple direction forward mode Taylor coefficient for op = ExpOp.
71
72The C++ source code corresponding to this operation is
73\verbatim
74        z = exp(x)
75\endverbatim
76
77\copydetails forward_unary1_op_dir
78*/
79template <class Base>
80inline void forward_exp_op_dir(
81        size_t q           ,
82        size_t r           ,
83        size_t i_z         ,
84        size_t i_x         ,
85        size_t cap_order   , 
86        Base*  taylor      )
87{       
88        // check assumptions
89        CPPAD_ASSERT_UNKNOWN( NumArg(ExpOp) == 1 );
90        CPPAD_ASSERT_UNKNOWN( NumRes(ExpOp) == 1 );
91        CPPAD_ASSERT_UNKNOWN( i_x < i_z );
92        CPPAD_ASSERT_UNKNOWN( q < cap_order );
93        CPPAD_ASSERT_UNKNOWN( 0 < q );
94
95        // Taylor coefficients corresponding to argument and result
96        size_t num_taylor_per_var = (cap_order-1) * r + 1;
97        Base* x = taylor + i_x * num_taylor_per_var;
98        Base* z = taylor + i_z * num_taylor_per_var; 
99
100        size_t m = (q-1)*r + 1;
101        for(size_t ell = 0; ell < r; ell++)
102        {       z[m+ell] = Base(q) * x[m+ell] * z[0];
103                for(size_t k = 1; k < q; k++)
104                        z[m+ell] += Base(k) * x[(k-1)*r+ell+1] * z[(q-k-1)*r+ell+1];
105                z[m+ell] /= Base(q);
106        }
107}
108
109/*!
110Zero order forward mode Taylor coefficient for result of op = ExpOp.
111
112The C++ source code corresponding to this operation is
113\verbatim
114        z = exp(x)
115\endverbatim
116
117\copydetails forward_unary1_op_0
118*/
119template <class Base>
120inline void forward_exp_op_0(
121        size_t i_z         ,
122        size_t i_x         ,
123        size_t cap_order   , 
124        Base*  taylor      )
125{
126        // check assumptions
127        CPPAD_ASSERT_UNKNOWN( NumArg(ExpOp) == 1 );
128        CPPAD_ASSERT_UNKNOWN( NumRes(ExpOp) == 1 );
129        CPPAD_ASSERT_UNKNOWN( i_x < i_z );
130        CPPAD_ASSERT_UNKNOWN( 0 < cap_order );
131
132        // Taylor coefficients corresponding to argument and result
133        Base* x = taylor + i_x * cap_order;
134        Base* z = taylor + i_z * cap_order;
135
136        z[0] = exp( x[0] );
137}
138/*!
139Reverse mode partial derivatives for result of op = ExpOp.
140
141The C++ source code corresponding to this operation is
142\verbatim
143        z = exp(x)
144\endverbatim
145
146\copydetails reverse_unary1_op
147*/
148
149template <class Base>
150inline void reverse_exp_op(
151        size_t      d            ,
152        size_t      i_z          ,
153        size_t      i_x          ,
154        size_t      cap_order    , 
155        const Base* taylor       ,
156        size_t      nc_partial   ,
157        Base*       partial      )
158{
159        // check assumptions
160        CPPAD_ASSERT_UNKNOWN( NumArg(ExpOp) == 1 );
161        CPPAD_ASSERT_UNKNOWN( NumRes(ExpOp) == 1 );
162        CPPAD_ASSERT_UNKNOWN( i_x < i_z );
163        CPPAD_ASSERT_UNKNOWN( d < cap_order );
164        CPPAD_ASSERT_UNKNOWN( d < nc_partial );
165
166        // Taylor coefficients and partials corresponding to argument
167        const Base* x  = taylor  + i_x * cap_order;
168        Base* px       = partial + i_x * nc_partial;
169
170        // Taylor coefficients and partials corresponding to result
171        const Base* z  = taylor  + i_z * cap_order;
172        Base* pz       = partial + i_z * nc_partial;
173
174        // loop through orders in reverse
175        size_t j, k;
176        j = d;
177        while(j)
178        {       // scale partial w.r.t z[j]
179                pz[j] /= Base(j);
180
181                for(k = 1; k <= j; k++)
182                {       px[k]   += pz[j] * Base(k) * z[j-k];   
183                        pz[j-k] += pz[j] * Base(k) * x[k];
184                }
185                --j;
186        }
187        px[0] += pz[0] * z[0];
188}
189
190} // END_CPPAD_NAMESPACE
191# endif
Note: See TracBrowser for help on using the repository browser.