CppAD: A C++ Algorithmic Differentiation Package  20130918
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-12 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Eclipse 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 $icode%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 $cref/AD function/glossary/AD Function/$$ corresponding to $icode 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 $icode i$$ and $icode j$$.
00053 
00054 $head f$$
00055 The object $icode f$$ has prototype
00056 $codei%
00057      ADFun<%Base%> %f%
00058 %$$
00059 Note that the $cref ADFun$$ object $icode f$$ is not $code const$$
00060 (see $cref/RevTwo Uses Forward/RevTwo/RevTwo Uses Forward/$$ below).
00061 
00062 $head x$$
00063 The argument $icode x$$ has prototype
00064 $codei%
00065      const %VectorBase% &%x%
00066 %$$
00067 (see $cref/VectorBase/RevTwo/VectorBase/$$ below)
00068 and its size 
00069 must be equal to $icode n$$, the dimension of the
00070 $cref/domain/seq_property/Domain/$$ space for $icode f$$.
00071 It specifies
00072 that point at which to evaluate the partial derivatives listed above.
00073 
00074 $head i$$
00075 The argument $icode i$$ has prototype
00076 $codei%
00077      const %VectorSize_t% &%i%
00078 %$$
00079 (see $cref/VectorSize_t/RevTwo/VectorSize_t/$$ below)
00080 We use $icode p$$ to denote the size of the vector $icode i$$.
00081 All of the indices in $icode i$$ 
00082 must be less than $icode m$$, the dimension of the
00083 $cref/range/seq_property/Range/$$ space for $icode f$$; i.e.,
00084 for $latex \ell = 0 , \ldots , p-1$$, $latex i[ \ell ]  < m$$.
00085 
00086 $head j$$
00087 The argument $icode j$$ has prototype
00088 $codei%
00089      const %VectorSize_t% &%j%
00090 %$$
00091 (see $cref/VectorSize_t/RevTwo/VectorSize_t/$$ below)
00092 and its size must be equal to $icode p$$,
00093 the size of the vector $icode i$$.
00094 All of the indices in $icode j$$ 
00095 must be less than $icode n$$; i.e.,
00096 for $latex \ell = 0 , \ldots , p-1$$, $latex j[ \ell ]  < n$$.
00097 
00098 $head ddw$$
00099 The result $icode ddw$$ has prototype
00100 $codei%
00101      %VectorBase% %ddw%
00102 %$$
00103 (see $cref/VectorBase/RevTwo/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 $icode VectorBase$$ must be a $cref SimpleVector$$ class with
00116 $cref/elements of type Base/SimpleVector/Elements of Specified Type/$$.
00117 The routine $cref CheckSimpleVector$$ will generate an error message
00118 if this is not the case.
00119 
00120 $head VectorSize_t$$
00121 The type $icode VectorSize_t$$ must be a $cref SimpleVector$$ class with
00122 $cref/elements of type size_t/SimpleVector/Elements of Specified Type/$$.
00123 The routine $cref CheckSimpleVector$$ will generate an error message
00124 if this is not the case.
00125 
00126 $head RevTwo Uses Forward$$
00127 After each call to $cref Forward$$,
00128 the object $icode f$$ contains the corresponding 
00129 $cref/Taylor coefficients/glossary/Taylor Coefficient/$$.
00130 After a call to $code RevTwo$$,
00131 the zero order Taylor coefficients correspond to
00132 $icode%f%.Forward(0, %x%)%$$
00133 and the other coefficients are unspecified.
00134 
00135 $head Examples$$
00136 $children%
00137      example/rev_two.cpp
00138 %$$
00139 The routine 
00140 $cref/RevTwo/rev_two.cpp/$$ is both an example and test.
00141 It returns $code true$$, if it succeeds and $code false$$ otherwise.
00142 
00143 $end
00144 -----------------------------------------------------------------------------
00145 */
00146 
00147 //  BEGIN CppAD namespace
00148 namespace CppAD {
00149 
00150 template <typename Base>
00151 template <typename VectorBase, typename VectorSize_t>
00152 VectorBase ADFun<Base>::RevTwo(
00153      const VectorBase   &x,
00154      const VectorSize_t &i,
00155      const VectorSize_t &j)
00156 {    size_t i1;
00157      size_t j1;
00158      size_t k;
00159      size_t l;
00160 
00161      size_t n = Domain();
00162      size_t m = Range();
00163      size_t p = i.size();
00164 
00165      // check VectorBase is Simple Vector class with Base elements
00166      CheckSimpleVector<Base, VectorBase>();
00167 
00168      // check VectorSize_t is Simple Vector class with size_t elements
00169      CheckSimpleVector<size_t, VectorSize_t>();
00170 
00171      CPPAD_ASSERT_KNOWN(
00172           x.size() == n,
00173           "RevTwo: Length of x not equal domain dimension for f."
00174      ); 
00175      CPPAD_ASSERT_KNOWN(
00176           i.size() == j.size(),
00177           "RevTwo: Lenght of the i and j vectors are not equal."
00178      );
00179      // point at which we are evaluating the second partials
00180      Forward(0, x);
00181 
00182      // dimension the return value
00183      VectorBase ddw(n * p);
00184 
00185      // direction vector in argument space
00186      VectorBase dx(n);
00187      for(j1 = 0; j1 < n; j1++)
00188           dx[j1] = Base(0);
00189 
00190      // direction vector in range space
00191      VectorBase w(m);
00192      for(i1 = 0; i1 < m; i1++)
00193           w[i1] = Base(0);
00194 
00195      // place to hold the results of a reverse calculation
00196      VectorBase r(n * 2);
00197 
00198      // check the indices in i and j
00199      for(l = 0; l < p; l++)
00200      {    i1 = i[l];
00201           j1 = j[l];
00202           CPPAD_ASSERT_KNOWN(
00203           i1 < m,
00204           "RevTwo: an eleemnt of i not less than range dimension for f."
00205           );
00206           CPPAD_ASSERT_KNOWN(
00207           j1 < n,
00208           "RevTwo: an element of j not less than domain dimension for f."
00209           );
00210      }
00211 
00212      // loop over all forward directions 
00213      for(j1 = 0; j1 < n; j1++) 
00214      {    // first order forward mode calculation done 
00215           bool first_done = false;
00216           for(l = 0; l < p; l++) if( j[l] == j1 )
00217           {    if( ! first_done )
00218                {    first_done = true;
00219 
00220                     // first order forward mode in j1 direction
00221                     dx[j1] = Base(1);
00222                     Forward(1, dx);
00223                     dx[j1] = Base(0);
00224                }
00225                // execute a reverse in this component direction
00226                i1    = i[l];
00227                w[i1] = Base(1);
00228                r     = Reverse(2, w);
00229                w[i1] = Base(0);
00230 
00231                // place the reverse result in return value
00232                for(k = 0; k < n; k++)
00233                     ddw[k * p + l] = r[k * 2 + 1];
00234           }
00235      }
00236      return ddw;
00237 }
00238 
00239 } // END CppAD namespace
00240 
00241 # endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines