CppAD: A C++ Algorithmic Differentiation Package  20130918
checkpoint.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_CHECKPOINT_INCLUDED
00003 # define CPPAD_CHECKPOINT_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 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 namespace CppAD { // BEGIN_CPPAD_NAMESPACE
00017 /*!
00018 \file checkpoint.hpp
00019 defining checkpoint functions.
00020 */
00021 
00022 /*
00023 $begin checkpoint$$
00024 $spell
00025      cppad.hpp
00026      CppAD
00027      checkpoint
00028      checkpointing
00029      algo
00030      afun
00031      const
00032 $$
00033 
00034 $section Checkpointing Functions$$
00035 $index function, checkpoint$$
00036 $index checkpoint, function$$
00037 
00038 $head Syntax$$
00039 $codei%checkpoint<%Base%> %afun%(%name%, %algo%, %ax%, %ay%)
00040 %afun%.option(%option_value%)
00041 %algo%(%ax%, %ay%)
00042 %afun%(%ax%, %ay%)
00043 checkpoint<%Base%>::clear()%$$
00044 
00045 $head Purpose$$
00046 You can reduce the size of the tape and memory required for AD by 
00047 checkpointing functions of the form $latex y = f(x)$$ where
00048 $latex f : B^n \rightarrow B^m$$.
00049 
00050 $head Method$$
00051 The $code checkpoint$$ class is derived from $code atomic_base$$
00052 and makes this easy.
00053 It implements all the $code atomic_base$$
00054 $cref/virtual functions/atomic_base/Virtual Functions/$$
00055 and hence its source code $code cppad/local/checkpoint.hpp$$
00056 provides an example implementation of $cref atomic_base$$.
00057 The difference is that $code checkpoint.hpp$$ uses AD 
00058 instead of user provided derivatives.
00059 
00060 $head constructor$$
00061 The constructor 
00062 $codei%
00063      checkpoint<%Base%> %afun%(%name%, %algo%, %ax%, %ay%)
00064 %$$
00065 cannot be called in $cref/parallel/ta_in_parallel/$$ mode.
00066 In addition, you cannot currently be recording 
00067 $codei%AD<%Base%>%$$ operations when the constructor is called.
00068 This class is implemented as a derived class of
00069 $cref/atomic_base/atomic_ctor/atomic_base/$$ and hence 
00070 some of its error message will refer to $code atomic_base$$.
00071 
00072 $head Base$$
00073 The type $icode Base$$ specifies the base type for AD operations.
00074 
00075 $head ADVector$$
00076 The type $icode ADVector$$ must be a
00077 $cref/simple vector class/SimpleVector/$$ with elements of type
00078 $codei%AD<%Base%>%$$.
00079 
00080 $head name$$
00081 This $icode checkpoint$$ constructor argument has prototype
00082 $codei%
00083      const char* %name%
00084 %$$
00085 It is the name used for error reporting.
00086 The suggested value for $icode name$$ is $icode afun$$; i.e.,
00087 the same name as used for the function.
00088 
00089 $head ax$$
00090 This argument has prototype
00091 $codei%
00092      const %ADVector%& %ax%
00093 %$$
00094 and size must be equal to $icode n$$.
00095 It specifies vector $latex x \in B^n$$ 
00096 at which an $codei%AD<%Base%>%$$ version of 
00097 $latex y = f(x)$$ is to be evaluated.
00098 
00099 $head ay$$
00100 This argument has prototype
00101 $codei%
00102      %ADVector%& %ay%
00103 %$$
00104 Its input size must be equal to $icode m$$ and does not change.
00105 The input values of its elements do not matter.
00106 Upon return, it is an $codei%AD<%Base%>%$$ version of 
00107 $latex y = f(x)$$.
00108 
00109 $head option$$
00110 The $code option$$ syntax can be used to set the type of sparsity
00111 pattern used by $icode afun$$.
00112 This is an $codei%atomic_base<%Base%>%$$ function and its documentation
00113 can be found at $cref atomic_option$$.
00114 
00115 $head algo$$
00116 The type of $icode algo$$ is arbitrary, except for the fact that
00117 the syntax
00118 $codei%
00119      %algo%(%ax%, %ay%)
00120 %$$ 
00121 must evaluate the function $latex y = f(x)$$ using
00122 $codei%AD<%Base%>%$$ operations.
00123 In addition, we assume that the 
00124 $cref/operation sequence/glossary/Operation/Sequence/$$
00125 does not depend on the value of $icode ax$$.
00126 
00127 $head afun$$
00128 Given $icode ax$$ it computes the corresponding value of $icode ay$$
00129 using the operation sequence corresponding to $icode algo$$. 
00130 If $codei%AD<%Base%>%$$ operations are being recorded,
00131 it enters the computation as single operation in the recording
00132 see $cref/start recording/Independent/Start Recording/$$.
00133 (Currently each use of $icode afun$$ actually corresponds to
00134 $icode%m%+%n%+2%$$ operations and creates $icode m$$ new variables, 
00135 but this is not part of the CppAD specifications and my change.)
00136 
00137 $head clear$$
00138 The $code atomic_base$$ class holds onto static work space in order to
00139 increase speed by avoiding system memory allocation calls.
00140 This call makes to work space $cref/available/ta_available/$$ to
00141 for other uses by the same thread.
00142 This should be called when you are done using the 
00143 user atomic functions for a specific value of $icode Base$$.
00144 
00145 $subhead Restriction$$
00146 The $code clear$$ routine cannot be called
00147 while in $cref/parallel/ta_in_parallel/$$ execution mode.
00148 
00149 $children%
00150      example/atomic/checkpoint.cpp
00151 %$$
00152 $head Example$$
00153 The file $cref checkpoint.cpp$$ contains an example and test
00154 of these operations.
00155 It returns true if it succeeds and false if it fails.
00156 
00157 $end
00158 */
00159 template <class Base>
00160 class checkpoint : public atomic_base<Base> {
00161 private:
00162      ADFun<Base> f_;
00163 public:
00164      /*!
00165      Constructor of a checkpoint object
00166 
00167      \param name [in]
00168      is the user's name for the AD version of this atomic operation.
00169 
00170      \param algo [in/out]
00171      user routine that compute AD function values
00172      (not const because state may change during evaluation).
00173 
00174      \param ax [in]
00175      argument value where algo operation sequence is taped.
00176 
00177      \param ay [out]
00178      function value at specified argument value.
00179      */
00180      template <class Algo, class ADVector>
00181      checkpoint(const char* name, 
00182           Algo& algo, const ADVector& ax, ADVector& ay)
00183      : atomic_base<Base>(name)
00184      {    CheckSimpleVector< CppAD::AD<Base> , ADVector>();
00185 
00186           // make a copy of ax because Independent modifies AD information
00187           ADVector x_tmp(ax);
00188           // delcare x_tmp as the independent variables
00189           Independent(x_tmp);
00190           // record mapping from x_tmp to ay
00191           algo(x_tmp, ay); 
00192           // create function f_ : x -> y
00193           f_.Dependent(ay);
00194           // suppress checking for nan in f_ results
00195           // (see optimize documentation for atomic functions)
00196           f_.check_for_nan(false);
00197           // now optimize (we expect to use this function many times).
00198           f_.optimize();
00199      }
00200      /*!
00201      Implement the user call to <tt>afun(ax, ay)</tt>.
00202      
00203      \tparam ADVector
00204      A simple vector class with elements of type <code>AD<Base></code>.
00205      
00206      \param id
00207      optional parameter which must be zero if present.
00208      
00209      \param ax
00210      is the argument vector for this call,
00211      <tt>ax.size()</tt> determines the number of arguments.
00212      
00213      \param ay
00214      is the result vector for this call,
00215      <tt>ay.size()</tt> determines the number of results.
00216      */
00217      template <class ADVector>
00218      void operator()(const ADVector& ax, ADVector& ay, size_t id = 0)
00219      {    CPPAD_ASSERT_KNOWN(
00220                id == 0,
00221                "checkpoint: id is non-zero in afun(ax, ay, id)"
00222           );
00223           this->atomic_base<Base>::operator()(ax, ay, id);
00224      }
00225      /*!
00226      Link from user_atomic to forward mode 
00227 
00228      \copydetails atomic_base::forward
00229      */
00230      virtual bool forward(
00231           size_t                    p ,
00232           size_t                    q ,
00233           const vector<bool>&      vx , 
00234                 vector<bool>&      vy , 
00235           const vector<Base>&      tx ,
00236                 vector<Base>&      ty )
00237      {
00238           CPPAD_ASSERT_UNKNOWN( f_.size_var() > 0 );
00239           CPPAD_ASSERT_UNKNOWN( tx.size() % (q+1) == 0 );
00240           CPPAD_ASSERT_UNKNOWN( ty.size() % (q+1) == 0 );
00241           size_t n = tx.size() / (q+1);
00242           size_t m = ty.size() / (q+1);
00243           bool ok  = true;    
00244           size_t i, j;
00245 
00246           // 2DO: test both forward and reverse vy information
00247           if( vx.size() > 0 )
00248           {    //Compute Jacobian sparsity pattern.
00249                vector< std::set<size_t> > s(m);
00250                if( n <= m )
00251                {    vector< std::set<size_t> > r(n);
00252                     for(j = 0; j < n; j++)
00253                          r[j].insert(j);
00254                     s = f_.ForSparseJac(n, r);
00255                }
00256                else
00257                {    vector< std::set<size_t> > r(m);
00258                     for(i = 0; i < m; i++)
00259                          r[i].insert(i);
00260                     s = f_.RevSparseJac(m, r);
00261                }
00262                std::set<size_t>::const_iterator itr;
00263                for(i = 0; i < m; i++)
00264                {    vy[i] = false;
00265                     for(itr = s[i].begin(); itr != s[i].end(); itr++)
00266                     {    j = *itr;
00267                          assert( j < n );
00268                          // y[i] depends on the value of x[j]
00269                          vy[i] |= vx[j];
00270                     }
00271                }
00272           }
00273           ty = f_.Forward(q, tx);
00274 
00275           // no longer need the Taylor coefficients in f_
00276           // (have to reconstruct them every time)
00277           f_.capacity_order(0);
00278           return ok;
00279      }
00280      /*!
00281      Link from user_atomic to reverse mode 
00282 
00283      \copydetails atomic_base::reverse
00284      */
00285      virtual bool reverse(
00286           size_t                    q  ,
00287           const vector<Base>&       tx ,
00288           const vector<Base>&       ty ,
00289                 vector<Base>&       px ,
00290           const vector<Base>&       py )
00291      {
00292           CPPAD_ASSERT_UNKNOWN( f_.size_var() > 0 );
00293           CPPAD_ASSERT_UNKNOWN( tx.size() % (q+1) == 0 );
00294           CPPAD_ASSERT_UNKNOWN( ty.size() % (q+1) == 0 );
00295           bool ok  = true;    
00296 
00297           // put proper forward mode coefficients in f_
00298 # ifdef NDEBUG
00299           f_.Forward(q, tx);
00300 # else
00301           size_t n = tx.size() / (q+1);
00302           size_t m = ty.size() / (q+1);
00303           CPPAD_ASSERT_UNKNOWN( px.size() == n * (q+1) );
00304           CPPAD_ASSERT_UNKNOWN( py.size() == m * (q+1) );
00305           size_t i, j, k;
00306           //
00307           vector<Base> check_ty = f_.Forward(q, tx);
00308           for(i = 0; i < m; i++)
00309           {    for(k = 0; k <= q; k++)
00310                {    j = i * (q+1) + k;
00311                     CPPAD_ASSERT_UNKNOWN( check_ty[j] == ty[j] );
00312                }
00313           }
00314 # endif
00315           // now can run reverse mode
00316           px = f_.Reverse(q+1, py);
00317 
00318           // no longer need the Taylor coefficients in f_
00319           // (have to reconstruct them every time)
00320           f_.capacity_order(0);
00321           return ok;
00322      }
00323      /*!
00324      Link from user_atomic to forward sparse Jacobian 
00325 
00326      \copydetails atomic_base::for_sparse_jac
00327      */
00328      virtual bool for_sparse_jac(
00329           size_t                                  q  ,
00330           const vector< std::set<size_t> >&       r  ,
00331                 vector< std::set<size_t> >&       s  )
00332      {
00333           bool ok = true;
00334           s = f_.ForSparseJac(q, r);
00335 
00336           // no longer need the forward mode sparsity pattern
00337           // (have to reconstruct them every time)
00338           f_.size_forward_set(0);
00339           
00340           return ok; 
00341      }
00342      /*!
00343      Link from user_atomic to forward sparse Jacobian 
00344 
00345      \copydetails atomic_base::for_sparse_jac
00346      */
00347      virtual bool for_sparse_jac(
00348           size_t                                  q  ,
00349           const vector<bool>&                     r  ,
00350                 vector<bool>&                     s  )
00351      {
00352           bool ok = true;
00353           s = f_.ForSparseJac(q, r);
00354 
00355           // no longer need the forward mode sparsity pattern
00356           // (have to reconstruct them every time)
00357           f_.size_forward_bool(0);
00358           
00359           return ok; 
00360      }
00361      /*!
00362      Link from user_atomic to forward sparse Jacobian 
00363 
00364      \copydetails atomic_base::rev_sparse_jac
00365      */
00366      virtual bool rev_sparse_jac(
00367           size_t                                  q  ,
00368           const vector< std::set<size_t> >&       rt ,
00369                 vector< std::set<size_t> >&       st )
00370      {
00371           bool ok  = true;
00372 
00373           // compute rt
00374           // 2DO: remove need for nz_compare all the time. It is only really
00375           // necessary when optimizer calls this member function.
00376           bool transpose = true;
00377           bool nz_compare = true;
00378           st = f_.RevSparseJac(q, rt, transpose, nz_compare);
00379 
00380           return ok; 
00381      }
00382      /*!
00383      Link from user_atomic to forward sparse Jacobian 
00384 
00385      \copydetails atomic_base::rev_sparse_jac
00386      */
00387      virtual bool rev_sparse_jac(
00388           size_t                                  q  ,
00389           const vector<bool>&                     rt ,
00390                 vector<bool>&                     st )
00391      {
00392           bool ok  = true;
00393 
00394           // compute rt
00395           bool transpose  = true;
00396           bool nz_compare = true;
00397           // 2DO: remove need for nz_compare all the time. It is only really
00398           // necessary when optimizer calls this member function.
00399           st = f_.RevSparseJac(q, rt, transpose, nz_compare);
00400 
00401           return ok; 
00402      }
00403      /*!
00404      Link from user_atomic to forward sparse Jacobian 
00405 
00406      \copydetails atomic_base::rev_sparse_hes
00407      */
00408      virtual bool rev_sparse_hes(
00409           const vector<bool>&                     vx ,
00410           const vector<bool>&                     s  ,
00411                 vector<bool>&                     t  ,
00412           size_t                                  q  ,
00413           const vector< std::set<size_t> >&       r  ,
00414           const vector< std::set<size_t> >&       u  ,
00415                 vector< std::set<size_t> >&       v  )
00416      {    size_t n       = v.size();
00417           size_t m       = u.size();
00418           CPPAD_ASSERT_UNKNOWN( r.size() == v.size() );
00419           CPPAD_ASSERT_UNKNOWN( s.size() == m );
00420           CPPAD_ASSERT_UNKNOWN( t.size() == n );
00421           bool ok        = true;
00422           bool transpose = true;
00423           std::set<size_t>::const_iterator itr;
00424 
00425           // compute sparsity pattern for T(x) = S(x) * f'(x)
00426           t = f_.RevSparseJac(1, s);
00427 # ifndef NDEBUG
00428           for(size_t j = 0; j < n; j++)
00429                CPPAD_ASSERT_UNKNOWN( vx[j] || ! t[j] )
00430 # endif
00431 
00432           // V(x) = f'(x)^T * g''(y) * f'(x) * R  +  g'(y) * f''(x) * R 
00433           // U(x) = g''(y) * f'(x) * R
00434           // S(x) = g'(y)
00435           
00436           // compute sparsity pattern for A(x) = f'(x)^T * U(x)
00437           vector< std::set<size_t> > a(n);
00438           a = f_.RevSparseJac(q, u, transpose);
00439 
00440           // set version of s
00441           vector< std::set<size_t> > set_s(1);
00442           CPPAD_ASSERT_UNKNOWN( set_s[0].empty() );
00443           size_t i;
00444           for(i = 0; i < m; i++)
00445                if( s[i] )
00446                     set_s[0].insert(i);
00447 
00448           // compute sparsity pattern for H(x) = (S(x) * F)''(x) * R
00449           // (store it in v)
00450           f_.ForSparseJac(q, r);
00451           v = f_.RevSparseHes(q, set_s, transpose);
00452 
00453           // compute sparsity pattern for V(x) = A(x) + H(x)
00454           for(i = 0; i < n; i++)
00455           {    for(itr = a[i].begin(); itr != a[i].end(); itr++)
00456                {    size_t j = *itr;
00457                     CPPAD_ASSERT_UNKNOWN( j < q );
00458                     v[i].insert(j);
00459                }
00460           }
00461 
00462           // no longer need the forward mode sparsity pattern
00463           // (have to reconstruct them every time)
00464           f_.size_forward_set(0);
00465 
00466           return ok;
00467      }
00468      /*!
00469      Link from user_atomic to forward sparse Jacobian 
00470 
00471      \copydetails atomic_base::rev_sparse_hes
00472      */
00473      virtual bool rev_sparse_hes(
00474           const vector<bool>&                     vx ,
00475           const vector<bool>&                     s  ,
00476                 vector<bool>&                     t  ,
00477           size_t                                  q  ,
00478           const vector<bool>&                     r  ,
00479           const vector<bool>&                     u  ,
00480                 vector<bool>&                     v  )
00481      {
00482           CPPAD_ASSERT_UNKNOWN( r.size() == v.size() );
00483           CPPAD_ASSERT_UNKNOWN( s.size() == u.size() / q );
00484           CPPAD_ASSERT_UNKNOWN( t.size() == v.size() / q );
00485           size_t n       = t.size();
00486           bool ok        = true;
00487           bool transpose = true;
00488           std::set<size_t>::const_iterator itr;
00489           size_t i, j;
00490 
00491           // compute sparsity pattern for T(x) = S(x) * f'(x)
00492           t = f_.RevSparseJac(1, s);
00493 # ifndef NDEBUG
00494           for(j = 0; j < n; j++)
00495                CPPAD_ASSERT_UNKNOWN( vx[j] || ! t[j] )
00496 # endif
00497 
00498           // V(x) = f'(x)^T * g''(y) * f'(x) * R  +  g'(y) * f''(x) * R 
00499           // U(x) = g''(y) * f'(x) * R
00500           // S(x) = g'(y)
00501 
00502           // compute sparsity pattern for A(x) = f'(x)^T * U(x)
00503           vector<bool> a(n * q);
00504           a = f_.RevSparseJac(q, u, transpose);
00505 
00506           // compute sparsity pattern for H(x) =(S(x) * F)''(x) * R
00507           // (store it in v)
00508           f_.ForSparseJac(q, r);
00509           v = f_.RevSparseHes(q, s, transpose);
00510 
00511           // compute sparsity pattern for V(x) = A(x) + H(x)
00512           for(i = 0; i < n; i++)
00513           {    for(j = 0; j < q; j++)
00514                     v[ i * q + j ] |= a[ i * q + j];
00515           }
00516 
00517           // no longer need the forward mode sparsity pattern
00518           // (have to reconstruct them every time)
00519           f_.size_forward_set(0);
00520 
00521           return ok;
00522      }
00523 };
00524 
00525 } // END_CPPAD_NAMESPACE
00526 # endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines