CppAD: A C++ Algorithmic Differentiation Package 20110419
rev_two.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_REV_TWO_INCLUDED
00003 # define CPPAD_REV_TWO_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Common Public License Version 1.0.
00011 
00012 A copy of this license is included in the COPYING file of this distribution.
00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00014 -------------------------------------------------------------------------- */
00015 
00016 /*
00017 $begin RevTwo$$
00018 $spell
00019         ddw
00020         typename
00021         Taylor
00022         const
00023 $$
00024 
00025 
00026 $index partial, second order driver$$
00027 $index second, order partial driver$$
00028 $index driver, second order partial$$
00029 
00030 $index easy, partial$$
00031 $index driver, easy partial$$
00032 $index partial, easy$$
00033 
00034 
00035 $section Reverse Mode Second Partial Derivative Driver$$
00036 
00037 $head Syntax$$
00038 $syntax%%ddw% = %f%.RevTwo(%x%, %i%, %j%)%$$
00039 
00040 
00041 $head Purpose$$
00042 We use $latex F : B^n \rightarrow B^m$$ to denote the
00043 $xref/glossary/AD Function/AD function/$$ corresponding to $italic f$$.
00044 The syntax above sets 
00045 $latex \[
00046         ddw [ k * p + \ell ]
00047         =
00048         \DD{ F_{i[ \ell ]} }{ x_{j[ \ell ]} }{ x_k } (x) 
00049 \] $$
00050 for $latex k = 0 , \ldots , n-1$$
00051 and $latex \ell = 0 , \ldots , p$$,
00052 where $latex p$$ is the size of the vectors $italic i$$ and $italic j$$.
00053 
00054 $head f$$
00055 The object $italic f$$ has prototype
00056 $syntax%
00057         ADFun<%Base%> %f%
00058 %$$
00059 Note that the $xref/ADFun/$$ object $italic f$$ is not $code const$$
00060 (see $xref/RevTwo/RevTwo Uses Forward/RevTwo Uses Forward/$$ below).
00061 
00062 $head x$$
00063 The argument $italic x$$ has prototype
00064 $syntax%
00065         const %VectorBase% &%x%
00066 %$$
00067 (see $xref/RevTwo/VectorBase/VectorBase/$$ below)
00068 and its size 
00069 must be equal to $italic n$$, the dimension of the
00070 $xref/seq_property/Domain/domain/$$ space for $italic f$$.
00071 It specifies
00072 that point at which to evaluate the partial derivatives listed above.
00073 
00074 $head i$$
00075 The argument $italic i$$ has prototype
00076 $syntax%
00077         const %VectorSize_t% &%i%
00078 %$$
00079 (see $xref/RevTwo/VectorSize_t/VectorSize_t/$$ below)
00080 We use $italic p$$ to denote the size of the vector $italic i$$.
00081 All of the indices in $italic i$$ 
00082 must be less than $italic m$$, the dimension of the
00083 $xref/seq_property/Range/range/$$ space for $italic f$$; i.e.,
00084 for $latex \ell = 0 , \ldots , p-1$$, $latex i[ \ell ]  < m$$.
00085 
00086 $head j$$
00087 The argument $italic j$$ has prototype
00088 $syntax%
00089         const %VectorSize_t% &%j%
00090 %$$
00091 (see $xref/RevTwo/VectorSize_t/VectorSize_t/$$ below)
00092 and its size must be equal to $italic p$$,
00093 the size of the vector $italic i$$.
00094 All of the indices in $italic j$$ 
00095 must be less than $italic n$$; i.e.,
00096 for $latex \ell = 0 , \ldots , p-1$$, $latex j[ \ell ]  < n$$.
00097 
00098 $head ddw$$
00099 The result $italic ddw$$ has prototype
00100 $syntax%
00101         %VectorBase% %ddw%
00102 %$$
00103 (see $xref/RevTwo/VectorBase/VectorBase/$$ below)
00104 and its size is $latex n * p$$.
00105 It contains the requested partial derivatives; to be specific,
00106 for $latex k = 0 , \ldots , n - 1 $$ 
00107 and $latex \ell = 0 , \ldots , p - 1$$
00108 $latex \[
00109         ddw [ k * p + \ell ]
00110         =
00111         \DD{ F_{i[ \ell ]} }{ x_{j[ \ell ]} }{ x_k } (x) 
00112 \] $$
00113 
00114 $head VectorBase$$
00115 The type $italic VectorBase$$ must be a $xref/SimpleVector/$$ class with
00116 $xref/SimpleVector/Elements of Specified Type/elements of type Base/$$.
00117 The routine $xref/CheckSimpleVector/$$ will generate an error message
00118 if this is not the case.
00119 
00120 $head VectorSize_t$$
00121 The type $italic VectorSize_t$$ must be a $xref/SimpleVector/$$ class with
00122 $xref/SimpleVector/Elements of Specified Type/elements of type size_t/$$.
00123 The routine $xref/CheckSimpleVector/$$ will generate an error message
00124 if this is not the case.
00125 
00126 $head RevTwo Uses Forward$$
00127 After each call to $xref/Forward/$$,
00128 the object $italic f$$ contains the corresponding 
00129 $xref/glossary/Taylor Coefficient/Taylor coefficients/$$.
00130 After $code RevTwo$$,
00131 the previous calls to $xref/Forward/$$ are undefined.
00132 
00133 $head Examples$$
00134 $children%
00135         example/rev_two.cpp
00136 %$$
00137 The routine 
00138 $xref/RevTwo.cpp//RevTwo/$$ is both an example and test.
00139 It returns $code true$$, if it succeeds and $code false$$ otherwise.
00140 
00141 $end
00142 -----------------------------------------------------------------------------
00143 */
00144 
00145 //  BEGIN CppAD namespace
00146 namespace CppAD {
00147 
00148 template <typename Base>
00149 template <typename VectorBase, typename VectorSize_t>
00150 VectorBase ADFun<Base>::RevTwo(
00151         const VectorBase   &x,
00152         const VectorSize_t &i,
00153         const VectorSize_t &j)
00154 {       size_t i1;
00155         size_t j1;
00156         size_t k;
00157         size_t l;
00158 
00159         size_t n = Domain();
00160         size_t m = Range();
00161         size_t p = i.size();
00162 
00163         // check VectorBase is Simple Vector class with Base elements
00164         CheckSimpleVector<Base, VectorBase>();
00165 
00166         // check VectorSize_t is Simple Vector class with size_t elements
00167         CheckSimpleVector<size_t, VectorSize_t>();
00168 
00169         CPPAD_ASSERT_KNOWN(
00170                 x.size() == n,
00171                 "RevTwo: Length of x not equal domain dimension for f."
00172         ); 
00173         CPPAD_ASSERT_KNOWN(
00174                 i.size() == j.size(),
00175                 "RevTwo: Lenght of the i and j vectors are not equal."
00176         );
00177         // point at which we are evaluating the second partials
00178         Forward(0, x);
00179 
00180         // dimension the return value
00181         VectorBase ddw(n * p);
00182 
00183         // direction vector in argument space
00184         VectorBase dx(n);
00185         for(j1 = 0; j1 < n; j1++)
00186                 dx[j1] = Base(0);
00187 
00188         // direction vector in range space
00189         VectorBase w(m);
00190         for(i1 = 0; i1 < m; i1++)
00191                 w[i1] = Base(0);
00192 
00193         // place to hold the results of a reverse calculation
00194         VectorBase r(n * 2);
00195 
00196         // check the indices in i and j
00197         for(l = 0; l < p; l++)
00198         {       i1 = i[l];
00199                 j1 = j[l];
00200                 CPPAD_ASSERT_KNOWN(
00201                 i1 < m,
00202                 "RevTwo: an eleemnt of i not less than range dimension for f."
00203                 );
00204                 CPPAD_ASSERT_KNOWN(
00205                 j1 < n,
00206                 "RevTwo: an element of j not less than domain dimension for f."
00207                 );
00208         }
00209 
00210         // loop over all forward directions 
00211         for(j1 = 0; j1 < n; j1++) 
00212         {       // first order forward mode calculation done 
00213                 bool first_done = false;
00214                 for(l = 0; l < p; l++) if( j[l] == j1 )
00215                 {       if( ! first_done )
00216                         {       first_done = true;
00217 
00218                                 // first order forward mode in j1 direction
00219                                 dx[j1] = Base(1);
00220                                 Forward(1, dx);
00221                                 dx[j1] = Base(0);
00222                         }
00223                         // execute a reverse in this component direction
00224                         i1    = i[l];
00225                         w[i1] = Base(1);
00226                         r     = Reverse(2, w);
00227                         w[i1] = Base(0);
00228 
00229                         // place the reverse result in return value
00230                         for(k = 0; k < n; k++)
00231                                 ddw[k * p + l] = r[k * 2 + 1];
00232                 }
00233         }
00234         return ddw;
00235 }
00236 
00237 } // END CppAD namespace
00238 
00239 # endif