subgraph_reverse.hpp

Go to the documentation of this file.

# ifndef CPPAD_CORE_SUBGRAPH_REVERSE_HPP
# define CPPAD_CORE_SUBGRAPH_REVERSE_HPP
/* --------------------------------------------------------------------------
CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell

CppAD is distributed under multiple licenses. This distribution is under
the terms of the
                    Eclipse Public License Version 1.0.

A copy of this license is included in the COPYING file of this distribution.
Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
-------------------------------------------------------------------------- */
/*
$begin subgraph_reverse$$
$spell
     resize
     subgraph
     Subgraphs
     dw
     Taylor
     Bool
     const
$$

$section Reverse Mode Using Subgraphs$$

$head Syntax$$
$icode%f%.subgraph_reverse(%select_domain%)
%$$
$icode%f%.subgraph_reverse(%q%, %ell%, %col%, %dw%)
%$$

$head Purpose$$
We use $latex F : B^n \rightarrow B^m$$ to denote the
$cref/AD function/glossary/AD Function/$$ corresponding to $icode f$$.
Reverse mode computes the derivative of the $cref Forward$$ mode
$cref/Taylor coefficients/glossary/Taylor Coefficient/$$
with respect to the domain variable $latex x$$.

$head Notation$$
We use the reverse mode
$cref/notation/reverse_any/Notation/$$ with the following change:
the vector
$cref/w^(k)/reverse_any/Notation/w^(k)/$$ is defined
$latex \[
w_i^{(k)} = \left\{ \begin{array}{ll}
     1 & {\rm if} \; k = q-1 \; \R{and} \; i = \ell
     \\
     0       & {\rm otherwise}
\end{array} \right.
\] $$

$head BaseVector$$
The type $icode BaseVector$$ must be a $cref SimpleVector$$ class with
$cref/elements of type/SimpleVector/Elements of Specified Type/$$
$icode Base$$.
The routine $cref CheckSimpleVector$$ will generate an error message
if this is not the case.

$head BoolVector$$
The type $icode BoolVector$$ is a $cref SimpleVector$$ class with
$cref/elements of type/SimpleVector/Elements of Specified Type/$$
$code bool$$.

$head SizeVector$$
The type $icode SizeVector$$ is a $cref SimpleVector$$ class with
$cref/elements of type/SimpleVector/Elements of Specified Type/$$
$code size_t$$.

$head select_domain$$
The argument $icode select_domain$$ has prototype
$codei%
     const %BoolVector%& %select_domain%
%$$
It has size $latex n$$ and specifies which independent variables
to include in the calculation.
If $icode%select_domain%[%j%]%$$ is false,
it is assumed that $latex u^{(k)}_j = 0$$ for $latex k > 0$$; i.e.,
the $th j$$ component of the Taylor coefficient for $latex x$$,
with order greater that zero, are zero; see
$cref/u^(k)/reverse_any/Notation/u^(k)/$$.

$head q$$
The argument $icode q$$ has prototype
$codei%
     size_t %q%
%$$
and specifies the number of Taylor coefficient orders to be differentiated.

$head ell$$
The argument $icode ell$$ has prototype
$codei%
     size_t %ell%
%$$
and specifies the dependent variable index that we are computing
the derivatives for; i.e. $latex \ell$$.
This index can only be used once per, and after, a call that selects
the independent variables using $icode select_domain$$.

$head col$$
This argument $icode col$$ has prototype
$codei%
     %SizeVector% %col%
%$$
The input size and value of its elements do not matter.
The $icode%col%.resize%$$ member function is used to change its size
to the number the number of possible non-zero derivative components.
For each $icode c$$,
$codei%
     %select_domain%[ %col%[%c%] ] == true
     %col%[%c%+1] >= %col%[%c%]
%$$
and the derivative with respect to the $th j$$ independent
variable is possibly non-zero where
$icode%j% = %col%[%c%]%$$.

$head dw$$
The argument $icode dw$$ has prototype
$codei%
     %Vector% %dw%
%$$
Its input size and value does not matter.
Upon return,
it is a vector with size $latex n \times q$$.
For $latex c = 0 , \ldots , %col%.size()-1$$,
and $latex k = 0, \ldots , q-1$$,
$latex \[
     dw[ j * q + k ] = W^{(1)} ( x )_{j,k}
\] $$
is the derivative of the specified Taylor coefficients w.r.t the $th j$$
independent variable where $icode%j% = %col%[%c%]%$$.
Note that this corresponds to the $cref reverse_any$$ convention when
$cref/w/reverse_any/w/$$ has size $icode%m% * %q%$$.

$head Example$$
$children%
     example/sparse/subgraph_reverse.cpp
%$$
The file
$cref subgraph_reverse.cpp$$
contains an example and test of this operation.
It returns true if it succeeds and false otherwise.

$end
*/
namespace CppAD { // BEGIN_CPPAD_NAMESPACE
/*!
\file subgraph_reverse.hpp
Compute derivatvies using reverse mode and subgraphs.
*/

/*!
Initialize reverse mode derivative computation on subgraphs.

\param select_domain
is a vector with size equal to the dimension of the domain for this function.
Only derivatives w.r.t. the components that are true will be computed.

\par subgraph_info_.map_user_op()
If the input size of this vector is zero,
its value for this player (play_) is computed.

\par subgraph_info.in_subgraph_
This vector is initialized for a reverse mode computation on subgraphs.

\par subgraph_info.select_domain()
This vector is set equal to the select_domain argument.

\par subgraph_info.process_range()
This vector is initialized to have size Range() and its elements are false.
*/

template <typename Base>
template <typename VectorBool>
void ADFun<Base>::subgraph_reverse( const VectorBool& select_domain )
{    using local::pod_vector;

     CPPAD_ASSERT_UNKNOWN(
          dep_taddr_.size() == subgraph_info_.n_dep()
     );
     CPPAD_ASSERT_UNKNOWN(
          size_t( select_domain.size() ) == subgraph_info_.n_ind()
     );

     // map_user_op
     if( subgraph_info_.map_user_op().size() == 0 )
          subgraph_info_.set_map_user_op(&play_);
     else
     {    CPPAD_ASSERT_UNKNOWN( subgraph_info_.check_map_user_op(&play_) );
     }
     CPPAD_ASSERT_UNKNOWN(
          subgraph_info_.map_user_op().size() == play_.num_op_rec()
     );

     // initialize for reverse mode subgraph computations
     subgraph_info_.init_rev(&play_, select_domain);
     CPPAD_ASSERT_UNKNOWN(
          subgraph_info_.in_subgraph().size() == play_.num_op_rec()
     );

     return;
}


/*!
Use reverse mode to compute derivative of Taylor coefficients on a subgraph.

The function
\f$ X : {\bf R} \times {\bf R}^{n \times q} \rightarrow {\bf R} \f$
is defined by
\f[
X(t , u) = \sum_{k=0}^{q-1} u^{(k)} t^k
\f]
The function
\f$ Y : {\bf R} \times {\bf R}^{n \times q} \rightarrow {\bf R} \f$
is defined by
\f[
Y(t , u) = F[ X(t, u) ]
\f]
The function
\f$ W : {\bf R}^{n \times q} \rightarrow {\bf R} \f$ is defined by
\f[
W(u) = \sum_{k=0}^{q-1} ( w^{(k)} )^{\rm T}
\frac{1}{k !} \frac{ \partial^k } { t^k } Y(0, u)
\f]

\param q
is the number of Taylor coefficient we are differentiating.

\param ell
is the component of the range that is selected for differentiation.

\param col
is the set of indices j = col[c] where the return value is defined.
If an index j is not in col, then either its derivative is zero,
or it is not in select_domain.

\param dw
Is a vector \f$ dw \f$ such that
for j = col[c],
\f$ k = 0 , \ldots , q-1 \f$
\f[
     dw[ j * q + k ] = W^{(1)} ( x )_{j,k}
\f]
where the matrix \f$ x \f$ is the value for \f$ u \f$
that corresponding to the forward mode Taylor coefficients
for the independent variables as specified by previous calls to Forward.

\par subgraph_info.process_range()
The element process_range[ell] is set to true by this operation.

\par subgraph_info.in_subgraph_
some of the elements of this vector are set to have value ell
(so it can not longer be used to determine the subgraph corresponding to
the ell-th dependent variable).
*/
template <typename Base>
template <typename VectorBase, typename SizeVector>
void ADFun<Base>::subgraph_reverse(
     size_t      q   ,
     size_t      ell ,
     SizeVector& col ,
     VectorBase& dw  )
{    using local::pod_vector;

     // check VectorBase is Simple Vector class with Base type elements
     CheckSimpleVector<Base, VectorBase>();
     CPPAD_ASSERT_KNOWN(
          q > 0,
          "The second argument to Reverse must be greater than zero."
     );
     CPPAD_ASSERT_KNOWN(
          num_order_taylor_ >= q,
          "Less than q Taylor coefficients are currently stored"
          " in this ADFun object."
     );
     CPPAD_ASSERT_KNOWN(
          num_direction_taylor_ == 1,
          "reverse mode for Forward(q, r, xq) with more than one direction"
          "\n(r > 1) is not yet supported."
     );
     CPPAD_ASSERT_KNOWN(
          ell < dep_taddr_.size(),
          "dependent variable index in to large for this function"
     );
     CPPAD_ASSERT_KNOWN(
          subgraph_info_.process_range()[ell] == false,
          "This dependent variable index has already been processed\n"
          "after the previous subgraph_reverse(select_domain)."
     );

     // subgraph of operators connected to dependent variable ell
     pod_vector<addr_t> subgraph;
     subgraph_info_.get_rev(
          &play_, dep_taddr_, addr_t(ell), subgraph
     );

     // Add all the atomic function call operators
     // for calls that have first operator in the subgraph
     local::subgraph::entire_call(&play_, subgraph);

     // sort the subgraph
     std::sort( subgraph.data(), subgraph.data() + subgraph.size() );

     /*
     // Use this printout for debugging
     std::cout << "{ ";
     for(size_t k = 0; k < subgraph.size(); k++)
     {    if( k > 0 )
               std::cout << ", ";
          std::cout << subgraph[k];
     }
     std::cout << "}\n";
     */

     // initialize subgraph_partial_ matrix to zero on subgraph
     Base zero(0);
     subgraph_partial_.resize(num_var_tape_ * q);
     for(size_t k = 0; k < subgraph.size(); ++k)
     {
          size_t               i_op = size_t( subgraph[k] );
          local::OpCode        op;
          const addr_t*        arg;
          size_t               i_var;
          play_.get_op_info(i_op, op, arg, i_var);
          if( NumRes(op) == 0 )
          {    CPPAD_ASSERT_UNKNOWN(
                    op == local::UserOp  ||
                    op == local::UsrapOp ||
                    op == local::UsravOp ||
                    op == local::UsrrpOp
               );
          }
          else
          {    CPPAD_ASSERT_UNKNOWN( i_var >= NumRes(op) );
               size_t j_var = i_var + 1 - NumRes(op);
               for(size_t i = j_var; i <= i_var; ++i)
               {    for(size_t j = 0; j < q; ++j)
                         subgraph_partial_[i * q + j] = zero;
               }
          }
     }

     // set partial to one for component we are differentiating
     subgraph_partial_[ dep_taddr_[ell] * q + q - 1] = Base(1);

     // evaluate the derivatives
     CPPAD_ASSERT_UNKNOWN( cskip_op_.size() == play_.num_op_rec() );
     CPPAD_ASSERT_UNKNOWN( load_op_.size()  == play_.num_load_op_rec() );
     size_t n = Domain();
     local::reverse_sweep(
          q - 1,
          n,
          num_var_tape_,
          &play_,
          cap_order_taylor_,
          taylor_.data(),
          q,
          subgraph_partial_.data(),
          cskip_op_.data(),
          load_op_,
          subgraph
     );

     // number of non-zero in return value
     size_t col_size       = 0;
     size_t subgraph_index = 0;
     while( subgraph_index < subgraph.size() )
     {    // check for InvOp
          if( subgraph[subgraph_index] > addr_t(n) )
               subgraph_index = subgraph.size();
          else
          {    ++col_size;
               ++subgraph_index;
          }
     }
     col.resize(col_size);

     // return the derivative values
     dw.resize(n * q);
     for(size_t c = 0; c < col_size; ++c)
     {    size_t i_op = subgraph[c];
          CPPAD_ASSERT_UNKNOWN( play_.GetOp(i_op) == local::InvOp );
          //
          size_t j = i_op - 1;
          CPPAD_ASSERT_UNKNOWN( i_op == play_.var2op(ind_taddr_[j]) );
          //
          // return paritial for this independent variable
          col[c] = j;
          for(size_t k = 0; k < q; k++)
               dw[j * q + k ] = subgraph_partial_[ind_taddr_[j] * q + k];
     }
     //
     CPPAD_ASSERT_KNOWN( ! ( hasnan(dw) && check_for_nan_ ) ,
          "f.subgraph_reverse(dw, q, ell): dw has a nan,\n"
          "but none of f's Taylor coefficents are nan."
     );
     //
     return;
}

} // END_CPPAD_NAMESPACE
# endif