CppAD: A C++ Algorithmic Differentiation Package  20130918
for_two.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_FOR_TWO_INCLUDED
00003 # define CPPAD_FOR_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 ForTwo$$
00018 $spell
00019      ddy
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 Forward Mode Second Partial Derivative Driver$$
00036 
00037 $head Syntax$$
00038 $icode%ddy% = %f%.ForTwo(%x%, %j%, %k%)%$$
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      ddy [ i * p + \ell ]
00047      =
00048      \DD{ F_i }{ x_{j[ \ell ]} }{ x_{k[ \ell ]} } (x) 
00049 \] $$
00050 for $latex i = 0 , \ldots , m-1$$
00051 and $latex \ell = 0 , \ldots , p$$,
00052 where $latex p$$ is the size of the vectors $icode j$$ and $icode k$$.
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/ForTwo Uses Forward/ForTwo/ForTwo 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/ForTwo/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 j$$
00075 The argument $icode j$$ has prototype
00076 $codei%
00077      const %VectorSize_t% &%j%
00078 %$$
00079 (see $cref/VectorSize_t/ForTwo/VectorSize_t/$$ below)
00080 We use $icode p$$ to denote the size of the vector $icode j$$.
00081 All of the indices in $icode j$$ 
00082 must be less than $icode n$$; i.e.,
00083 for $latex \ell = 0 , \ldots , p-1$$, $latex j[ \ell ]  < n$$.
00084 
00085 $head k$$
00086 The argument $icode k$$ has prototype
00087 $codei%
00088      const %VectorSize_t% &%k%
00089 %$$
00090 (see $cref/VectorSize_t/ForTwo/VectorSize_t/$$ below)
00091 and its size must be equal to $icode p$$,
00092 the size of the vector $icode j$$.
00093 All of the indices in $icode k$$ 
00094 must be less than $icode n$$; i.e.,
00095 for $latex \ell = 0 , \ldots , p-1$$, $latex k[ \ell ]  < n$$.
00096 
00097 $head ddy$$
00098 The result $icode ddy$$ has prototype
00099 $codei%
00100      %VectorBase% %ddy%
00101 %$$
00102 (see $cref/VectorBase/ForTwo/VectorBase/$$ below)
00103 and its size is $latex m * p$$.
00104 It contains the requested partial derivatives; to be specific,
00105 for $latex i = 0 , \ldots , m - 1 $$ 
00106 and $latex \ell = 0 , \ldots , p - 1$$
00107 $latex \[
00108      ddy [ i * p + \ell ]
00109      =
00110      \DD{ F_i }{ x_{j[ \ell ]} }{ x_{k[ \ell ]} } (x) 
00111 \] $$
00112 
00113 $head VectorBase$$
00114 The type $icode VectorBase$$ must be a $cref SimpleVector$$ class with
00115 $cref/elements of type Base/SimpleVector/Elements of Specified Type/$$.
00116 The routine $cref CheckSimpleVector$$ will generate an error message
00117 if this is not the case.
00118 
00119 $head VectorSize_t$$
00120 The type $icode VectorSize_t$$ must be a $cref SimpleVector$$ class with
00121 $cref/elements of type size_t/SimpleVector/Elements of Specified Type/$$.
00122 The routine $cref CheckSimpleVector$$ will generate an error message
00123 if this is not the case.
00124 
00125 $head ForTwo Uses Forward$$
00126 After each call to $cref Forward$$,
00127 the object $icode f$$ contains the corresponding 
00128 $cref/Taylor coefficients/glossary/Taylor Coefficient/$$.
00129 After a call to $code ForTwo$$,
00130 the zero order Taylor coefficients correspond to
00131 $icode%f%.Forward(0, %x%)%$$
00132 and the other coefficients are unspecified.
00133 
00134 $head Examples$$
00135 $children%
00136      example/for_two.cpp
00137 %$$
00138 The routine 
00139 $cref/ForTwo/for_two.cpp/$$ is both an example and test.
00140 It returns $code true$$, if it succeeds and $code false$$ otherwise.
00141 
00142 $end
00143 -----------------------------------------------------------------------------
00144 */
00145 
00146 //  BEGIN CppAD namespace
00147 namespace CppAD {
00148 
00149 template <typename Base>
00150 template <typename VectorBase, typename VectorSize_t>
00151 VectorBase ADFun<Base>::ForTwo(
00152      const VectorBase   &x, 
00153      const VectorSize_t &j,
00154      const VectorSize_t &k)
00155 {    size_t i;
00156      size_t j1;
00157      size_t k1;
00158      size_t l;
00159 
00160      size_t n = Domain();
00161      size_t m = Range();
00162      size_t p = j.size();
00163 
00164      // check VectorBase is Simple Vector class with Base type elements
00165      CheckSimpleVector<Base, VectorBase>();
00166 
00167      // check VectorSize_t is Simple Vector class with size_t elements
00168      CheckSimpleVector<size_t, VectorSize_t>();
00169 
00170      CPPAD_ASSERT_KNOWN(
00171           x.size() == n,
00172           "ForTwo: Length of x not equal domain dimension for f."
00173      ); 
00174      CPPAD_ASSERT_KNOWN(
00175           j.size() == k.size(),
00176           "ForTwo: Lenght of the j and k vectors are not equal."
00177      );
00178      // point at which we are evaluating the second partials
00179      Forward(0, x);
00180 
00181 
00182      // dimension the return value
00183      VectorBase ddy(m * p);
00184 
00185      // allocate memory to hold all possible diagonal Taylor coefficients
00186      // (for large sparse cases, this is not efficient)
00187      VectorBase D(m * n);
00188 
00189      // boolean flag for which diagonal coefficients are computed
00190      CppAD::vector<bool> c(n);
00191      for(j1 = 0; j1 < n; j1++)
00192           c[j1] = false;
00193 
00194      // direction vector in argument space
00195      VectorBase dx(n);
00196      for(j1 = 0; j1 < n; j1++)
00197           dx[j1] = Base(0);
00198 
00199      // result vector in range space
00200      VectorBase dy(m);
00201 
00202      // compute the diagonal coefficients that are needed
00203      for(l = 0; l < p; l++)
00204      {    j1 = j[l];
00205           k1 = k[l];
00206           CPPAD_ASSERT_KNOWN(
00207           j1 < n,
00208           "ForTwo: an element of j not less than domain dimension for f."
00209           );
00210           CPPAD_ASSERT_KNOWN(
00211           k1 < n,
00212           "ForTwo: an element of k not less than domain dimension for f."
00213           );
00214           size_t count = 2;
00215           while(count)
00216           {    count--;
00217                if( ! c[j1] )
00218                {    // diagonal term in j1 direction
00219                     c[j1]  = true;
00220                     dx[j1] = Base(1);
00221                     Forward(1, dx);
00222 
00223                     dx[j1] = Base(0);
00224                     dy     = Forward(2, dx);
00225                     for(i = 0; i < m; i++)
00226                          D[i * n + j1 ] = dy[i];
00227                } 
00228                j1 = k1;
00229           }
00230      }
00231      // compute all the requested cross partials
00232      for(l = 0; l < p; l++)
00233      {    j1 = j[l];
00234           k1 = k[l];
00235           if( j1 == k1 )
00236           {    for(i = 0; i < m; i++)
00237                     ddy[i * p + l] = Base(2) * D[i * n + j1];
00238           }
00239           else
00240           {
00241                // cross term in j1 and k1 directions
00242                dx[j1] = Base(1);
00243                dx[k1] = Base(1);
00244                Forward(1, dx);
00245 
00246                dx[j1] = Base(0);
00247                dx[k1] = Base(0);
00248                dy = Forward(2, dx);
00249 
00250                // place result in return value
00251                for(i = 0; i < m; i++)
00252                     ddy[i * p + l] = dy[i] - D[i*n+j1] - D[i*n+k1];
00253 
00254           }
00255      }
00256      return ddy;
00257 }
00258 
00259 } // END CppAD namespace
00260 
00261 # endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines