source: trunk/cppad/local/log_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.8 KB
Line 
1/* $Id: log_op.hpp 3301 2014-05-24 05:20:21Z bradbell $ */
2# ifndef CPPAD_LOG_OP_INCLUDED
3# define CPPAD_LOG_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
16namespace CppAD { // BEGIN_CPPAD_NAMESPACE
17/*!
18\file log_op.hpp
19Forward and reverse mode calculations for z = log(x).
20*/
21
22/*!
23Compute forward mode Taylor coefficient for result of op = LogOp.
24
25The C++ source code corresponding to this operation is
26\verbatim
27        z = log(x)
28\endverbatim
29
30\copydetails forward_unary1_op
31*/
32template <class Base>
33inline void forward_log_op(
34        size_t p           ,
35        size_t q           ,
36        size_t i_z         ,
37        size_t i_x         ,
38        size_t cap_order   , 
39        Base*  taylor      )
40{       
41        size_t k;
42
43        // check assumptions
44        CPPAD_ASSERT_UNKNOWN( NumArg(LogOp) == 1 );
45        CPPAD_ASSERT_UNKNOWN( NumRes(LogOp) == 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        if( p == 0 )
55        {       z[0] = log( x[0] );
56                p++;
57                if( q == 0 )
58                        return;
59        }
60        if ( p == 1 )
61        {       z[1] = x[1] / x[0];
62                p++;
63        }
64        for(size_t j = p; j <= q; j++)
65        {
66                z[j] = -z[1] * x[j-1];
67                for(k = 2; k < j; k++)
68                        z[j] -= Base(k) * z[k] * x[j-k];
69                z[j] /= Base(j);
70                z[j] += x[j];
71                z[j] /= x[0];
72        }
73}
74
75/*!
76Muiltiple directions Taylor coefficient for op = LogOp.
77
78The C++ source code corresponding to this operation is
79\verbatim
80        z = log(x)
81\endverbatim
82
83\copydetails forward_unary1_op_dir
84*/
85template <class Base>
86inline void forward_log_op_dir(
87        size_t q           ,
88        size_t r           ,
89        size_t i_z         ,
90        size_t i_x         ,
91        size_t cap_order   , 
92        Base*  taylor      )
93{       
94
95        // check assumptions
96        CPPAD_ASSERT_UNKNOWN( NumArg(LogOp) == 1 );
97        CPPAD_ASSERT_UNKNOWN( NumRes(LogOp) == 1 );
98        CPPAD_ASSERT_UNKNOWN( i_x < i_z );
99        CPPAD_ASSERT_UNKNOWN( 0 < q );
100        CPPAD_ASSERT_UNKNOWN( q < cap_order );
101
102        // Taylor coefficients corresponding to argument and result
103        size_t num_taylor_per_var = (cap_order-1) * r + 1;
104        Base* x = taylor + i_x * num_taylor_per_var;
105        Base* z = taylor + i_z * num_taylor_per_var;
106
107        size_t m = (q-1) * r + 1;
108        for(size_t ell = 0; ell < r; ell++)
109        {       z[m+ell] = Base(q) * x[m+ell];
110                for(size_t k = 1; k < q; k++)
111                        z[m+ell] -= Base(k) * z[(k-1)*r+1+ell] * x[(q-k-1)*r+1+ell];
112                z[m+ell] /= (Base(q) * x[0]);
113        }       
114}
115
116/*!
117Compute zero order forward mode Taylor coefficient for result of op = LogOp.
118
119The C++ source code corresponding to this operation is
120\verbatim
121        z = log(x)
122\endverbatim
123
124\copydetails forward_unary1_op_0
125*/
126template <class Base>
127inline void forward_log_op_0(
128        size_t i_z         ,
129        size_t i_x         ,
130        size_t cap_order   , 
131        Base*  taylor      )
132{
133
134        // check assumptions
135        CPPAD_ASSERT_UNKNOWN( NumArg(LogOp) == 1 );
136        CPPAD_ASSERT_UNKNOWN( NumRes(LogOp) == 1 );
137        CPPAD_ASSERT_UNKNOWN( i_x < i_z );
138        CPPAD_ASSERT_UNKNOWN( 0 < cap_order );
139
140        // Taylor coefficients corresponding to argument and result
141        Base* x = taylor + i_x * cap_order;
142        Base* z = taylor + i_z * cap_order;
143
144        z[0] = log( x[0] );
145}
146
147/*!
148Compute reverse mode partial derivatives for result of op = LogOp.
149
150The C++ source code corresponding to this operation is
151\verbatim
152        z = log(x)
153\endverbatim
154
155\copydetails reverse_unary1_op
156*/
157
158template <class Base>
159inline void reverse_log_op(
160        size_t      d            ,
161        size_t      i_z          ,
162        size_t      i_x          ,
163        size_t      cap_order    , 
164        const Base* taylor       ,
165        size_t      nc_partial   ,
166        Base*       partial      )
167{       size_t j, k;   
168
169        // check assumptions
170        CPPAD_ASSERT_UNKNOWN( NumArg(LogOp) == 1 );
171        CPPAD_ASSERT_UNKNOWN( NumRes(LogOp) == 1 );
172        CPPAD_ASSERT_UNKNOWN( i_x < i_z );
173        CPPAD_ASSERT_UNKNOWN( d < cap_order );
174        CPPAD_ASSERT_UNKNOWN( d < nc_partial );
175
176        // Taylor coefficients and partials corresponding to argument
177        const Base* x  = taylor  + i_x * cap_order;
178        Base* px       = partial + i_x * nc_partial;
179
180        // Taylor coefficients and partials corresponding to result
181        const Base* z  = taylor  + i_z * cap_order;
182        Base* pz       = partial + i_z * nc_partial;
183
184        j = d;
185        while(j)
186        {       // scale partial w.r.t z[j]
187                pz[j]   /= x[0];
188
189                px[0]   -= pz[j] * z[j];
190                px[j]   += pz[j];
191
192                // further scale partial w.r.t. z[j]
193                pz[j]   /= Base(j);
194
195                for(k = 1; k < j; k++)
196                {       pz[k]   -= pz[j] * Base(k) * x[j-k];
197                        px[j-k] -= pz[j] * Base(k) * z[k];
198                }
199                --j;
200        }
201        px[0] += pz[0] / x[0];
202}
203
204} // END_CPPAD_NAMESPACE
205# endif
Note: See TracBrowser for help on using the repository browser.