// $Id$ # ifndef CPPAD_CSUM_OP_INCLUDED # define CPPAD_CSUM_OP_INCLUDED /* -------------------------------------------------------------------------- CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 Bradley M. Bell CppAD is distributed under multiple licenses. This distribution is under the terms of the GNU General Public License Version 3. 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. -------------------------------------------------------------------------- */ namespace CppAD { // BEGIN_CPPAD_NAMESPACE /*! \file csum_op.hpp Forward, reverse and sparsity calculations for cummulative summation. */ /*! Compute forward mode Taylor coefficients for result of op = CsumOp. This operation is \verbatim z = s + x(1) + ... + x(m) - y(1) - ... - y(n). \endverbatim \tparam Base base type for the operator; i.e., this operation was recorded using AD< \a Base > and computations by this routine are done using type \a Base. \param p lowest order of the Taylor coefficient that we are computing. \param q highest order of the Taylor coefficient that we are computing. \param i_z variable index corresponding to the result for this operation; i.e. the row index in \a taylor corresponding to z. \param arg \a arg[0] is the number of addition variables in this cummulative summation; i.e., m. \n \a arg[1] is the number of subtraction variables in this cummulative summation; i.e., \c m. \n parameter[ arg[2] ] is the parameter value \c s in this cummunative summation. \n arg[2+i] for i = 1 , ... , m is the variable index of x(i). \n arg[2+arg[0]+i] for i = 1 , ... , n is the variable index of y(i). \param num_par is the number of parameters in \a parameter. \param parameter is the parameter vector for this operation sequence. \param cap_order number of colums in the matrix containing all the Taylor coefficients. \param taylor \b Input: taylor [ arg[2+i] * cap_order + k ] for i = 1 , ... , m and k = 0 , ... , q is the k-th order Taylor coefficient corresponding to x(i) \n \b Input: taylor [ arg[2+m+i] * cap_order + k ] for i = 1 , ... , n and k = 0 , ... , q is the k-th order Taylor coefficient corresponding to y(i) \n \b Input: taylor [ i_z * cap_order + k ] for k = 0 , ... , p, is the k-th order Taylor coefficient corresponding to z. \n \b Output: taylor [ i_z * cap_order + k ] for k = p , ... , q, is the \a k-th order Taylor coefficient corresponding to z. */ template inline void forward_csum_op( size_t p , size_t q , size_t i_z , const addr_t* arg , size_t num_par , const Base* parameter , size_t cap_order , Base* taylor ) { Base zero(0); size_t i, j, k; // check assumptions CPPAD_ASSERT_UNKNOWN( NumRes(CSumOp) == 1 ); CPPAD_ASSERT_UNKNOWN( q < cap_order ); CPPAD_ASSERT_UNKNOWN( p <= q ); CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < num_par ); CPPAD_ASSERT_UNKNOWN( arg[0] + arg[1] == arg[ arg[0] + arg[1] + 3 ] ); // Taylor coefficients corresponding to result Base* z = taylor + i_z * cap_order; for(k = p; k <= q; k++) z[k] = zero; if( p == 0 ) z[p] = parameter[ arg[2] ]; Base* x; i = arg[0]; j = 2; while(i--) { CPPAD_ASSERT_UNKNOWN( size_t(arg[j+1]) < i_z ); x = taylor + arg[++j] * cap_order; for(k = p; k <= q; k++) z[k] += x[k]; } i = arg[1]; while(i--) { CPPAD_ASSERT_UNKNOWN( size_t(arg[j+1]) < i_z ); x = taylor + arg[++j] * cap_order; for(k = p; k <= q; k++) z[k] -= x[k]; } } /*! Multiple direction forward mode Taylor coefficients for op = CsumOp. This operation is \verbatim z = s + x(1) + ... + x(m) - y(1) - ... - y(n). \endverbatim \tparam Base base type for the operator; i.e., this operation was recorded using AD and computations by this routine are done using type \a Base. \param q order ot the Taylor coefficients that we are computing. \param r number of directions for Taylor coefficients that we are computing. \param i_z variable index corresponding to the result for this operation; i.e. the row index in \a taylor corresponding to z. \param arg \a arg[0] is the number of addition variables in this cummulative summation; i.e., m. \n \a arg[1] is the number of subtraction variables in this cummulative summation; i.e., \c m. \n parameter[ arg[2] ] is the parameter value \c s in this cummunative summation. \n arg[2+i] for i = 1 , ... , m is the variable index of x(i). \n arg[2+arg[0]+i] for i = 1 , ... , n is the variable index of y(i). \param num_par is the number of parameters in \a parameter. \param parameter is the parameter vector for this operation sequence. \param cap_order number of colums in the matrix containing all the Taylor coefficients. \param taylor \b Input: taylor [ arg[2+i]*((cap_order-1)*r + 1) + 0 ] for i = 1 , ... , m is the 0-th order Taylor coefficient corresponding to x(i) and taylor [ arg[2+i]*((cap_order-1)*r + 1) + (q-1)*r + ell + 1 ] for i = 1 , ... , m, ell = 0 , ... , r-1 is the q-th order Taylor coefficient corresponding to x(i) and direction ell. \n \b Input: taylor [ arg[2+m+i]*((cap_order-1)*r + 1) + 0 ] for i = 1 , ... , n is the 0-th order Taylor coefficient corresponding to y(i) and taylor [ arg[2+m+i]*((cap_order-1)*r + 1) + (q-1)*r + ell + 1 ] for i = 1 , ... , n, ell = 0 , ... , r-1 is the q-th order Taylor coefficient corresponding to y(i) and direction ell. \n \b Output: taylor [ i_z*((cap_order-1)*r+1) + (q-1)*r + ell + 1 ] is the \a q-th order Taylor coefficient corresponding to z for direction ell = 0 , ... , r-1. */ template inline void forward_csum_op_dir( size_t q , size_t r , size_t i_z , const addr_t* arg , size_t num_par , const Base* parameter , size_t cap_order , Base* taylor ) { Base zero(0); size_t i, j, ell; // check assumptions CPPAD_ASSERT_UNKNOWN( NumRes(CSumOp) == 1 ); CPPAD_ASSERT_UNKNOWN( q < cap_order ); CPPAD_ASSERT_UNKNOWN( 0 < q ); CPPAD_ASSERT_UNKNOWN( size_t(arg[2]) < num_par ); CPPAD_ASSERT_UNKNOWN( arg[0] + arg[1] == arg[ arg[0] + arg[1] + 3 ] ); // Taylor coefficients corresponding to result size_t num_taylor_per_var = (cap_order-1) * r + 1; size_t m = (q-1)*r + 1; Base* z = taylor + i_z * num_taylor_per_var + m; for(ell = 0; ell < r; ell++) z[ell] = zero; Base* x; i = arg[0]; j = 2; while(i--) { CPPAD_ASSERT_UNKNOWN( size_t(arg[j+1]) < i_z ); x = taylor + arg[++j] * num_taylor_per_var + m; for(ell = 0; ell < r; ell++) z[ell] += x[ell]; } i = arg[1]; while(i--) { CPPAD_ASSERT_UNKNOWN( size_t(arg[j+1]) < i_z ); x = taylor + arg[++j] * num_taylor_per_var + m; for(ell = 0; ell < r; ell++) z[ell] -= x[ell]; } } /*! Compute reverse mode Taylor coefficients for result of op = CsumOp. This operation is \verbatim z = q + x(1) + ... + x(m) - y(1) - ... - y(n). H(y, x, w, ...) = G[ z(x, y), y, x, w, ... ] \endverbatim \tparam Base base type for the operator; i.e., this operation was recorded using AD< \a Base > and computations by this routine are done using type \a Base. \param d order the highest order Taylor coefficient that we are computing the partial derivatives with respect to. \param i_z variable index corresponding to the result for this operation; i.e. the row index in \a taylor corresponding to z. \param arg \a arg[0] is the number of addition variables in this cummulative summation; i.e., m. \n \a arg[1] is the number of subtraction variables in this cummulative summation; i.e., \c m. \n parameter[ arg[2] ] is the parameter value \c q in this cummunative summation. \n arg[2+i] for i = 1 , ... , m is the value x(i). \n arg[2+arg[0]+i] for i = 1 , ... , n is the value y(i). \param nc_partial number of colums in the matrix containing all the partial derivatives. \param partial \b Input: partial [ arg[2+i] * nc_partial + k ] for i = 1 , ... , m and k = 0 , ... , d is the partial derivative of G(z, y, x, w, ...) with respect to the k-th order Taylor coefficient corresponding to x(i) \n \b Input: partial [ arg[2+m+i] * nc_partial + k ] for i = 1 , ... , n and k = 0 , ... , d is the partial derivative of G(z, y, x, w, ...) with respect to the k-th order Taylor coefficient corresponding to y(i) \n \b Input: partial [ i_z * nc_partial + k ] for i = 1 , ... , n and k = 0 , ... , d is the partial derivative of G(z, y, x, w, ...) with respect to the k-th order Taylor coefficient corresponding to \c z. \n \b Output: partial [ arg[2+i] * nc_partial + k ] for i = 1 , ... , m and k = 0 , ... , d is the partial derivative of H(y, x, w, ...) with respect to the k-th order Taylor coefficient corresponding to x(i) \n \b Output: partial [ arg[2+m+i] * nc_partial + k ] for i = 1 , ... , n and k = 0 , ... , d is the partial derivative of H(y, x, w, ...) with respect to the k-th order Taylor coefficient corresponding to y(i) */ template inline void reverse_csum_op( size_t d , size_t i_z , const addr_t* arg , size_t nc_partial , Base* partial ) { // check assumptions CPPAD_ASSERT_UNKNOWN( NumRes(CSumOp) == 1 ); CPPAD_ASSERT_UNKNOWN( d < nc_partial ); // Taylor coefficients and partial derivative corresponding to result Base* pz = partial + i_z * nc_partial; Base* px; size_t i, j, k; size_t d1 = d + 1; i = arg[0]; j = 2; while(i--) { CPPAD_ASSERT_UNKNOWN( size_t(arg[j+1]) < i_z ); px = partial + arg[++j] * nc_partial; k = d1; while(k--) px[k] += pz[k]; } i = arg[1]; while(i--) { CPPAD_ASSERT_UNKNOWN( size_t(arg[j+1]) < i_z ); px = partial + arg[++j] * nc_partial; k = d1; while(k--) px[k] -= pz[k]; } } /*! Forward mode Jacobian sparsity pattern for CSumOp operator. This operation is \verbatim z = q + x(1) + ... + x(m) - y(1) - ... - y(n). \endverbatim \tparam Vector_set is the type used for vectors of sets. It can be either \c sparse_pack, \c sparse_set, or \c sparse_list. \param i_z variable index corresponding to the result for this operation; i.e. the index in \a sparsity corresponding to z. \param arg \a arg[0] is the number of addition variables in this cummulative summation; i.e., m + n. \n \a arg[1] is the number of subtraction variables in this cummulative summation; i.e., \c m. \n parameter[ arg[2] ] is the parameter value \c q in this cummunative summation. \n arg[2+i] for i = 1 , ... , m is the value x(i). \n arg[2+arg[1]+i] for i = 1 , ... , n is the value y(i). \param sparsity \b Input: For i = 1 , ... , m, the set with index \a arg[2+i] in \a sparsity is the sparsity bit pattern for x(i). This identifies which of the independent variables the variable x(i) depends on. \n \b Input: For i = 1 , ... , n, the set with index \a arg[2+arg[0]+i] in \a sparsity is the sparsity bit pattern for x(i). This identifies which of the independent variables the variable y(i) depends on. \n \b Output: The set with index \a i_z in \a sparsity is the sparsity bit pattern for z. This identifies which of the independent variables the variable z depends on. */ template inline void forward_sparse_jacobian_csum_op( size_t i_z , const addr_t* arg , Vector_set& sparsity ) { sparsity.clear(i_z); size_t i, j; i = arg[0] + arg[1]; j = 2; while(i--) { CPPAD_ASSERT_UNKNOWN( size_t(arg[j+1]) < i_z ); sparsity.binary_union( i_z , // index in sparsity for result i_z , // index in sparsity for left operand arg[++j] , // index for right operand sparsity // sparsity vector for right operand ); } } /*! Reverse mode Jacobian sparsity pattern for CSumOp operator. This operation is \verbatim z = q + x(1) + ... + x(m) - y(1) - ... - y(n). H(y, x, w, ...) = G[ z(x, y), y, x, w, ... ] \endverbatim \tparam Vector_set is the type used for vectors of sets. It can be either \c sparse_pack, \c sparse_set, or \c sparse_list. \param i_z variable index corresponding to the result for this operation; i.e. the index in \a sparsity corresponding to z. \param arg \a arg[0] is the number of addition variables in this cummulative summation; i.e., m + n. \n \a arg[1] is the number of subtraction variables in this cummulative summation; i.e., \c m. \n parameter[ arg[2] ] is the parameter value \c q in this cummunative summation. \n arg[2+i] for i = 1 , ... , m is the value x(i). \n arg[2+arg[1]+i] for i = 1 , ... , n is the value y(i). \param sparsity For i = 1 , ... , m, the set with index \a arg[2+i] in \a sparsity is the sparsity bit pattern for x(i). This identifies which of the dependent variables depend on x(i). On input, the sparsity patter corresponds to \c G, and on ouput it corresponds to \c H. \n For i = 1 , ... , m, the set with index \a arg[2+arg[0]+i] in \a sparsity is the sparsity bit pattern for y(i). This identifies which of the dependent variables depend on y(i). On input, the sparsity patter corresponds to \c G, and on ouput it corresponds to \c H. \n \b Input: The set with index \a i_z in \a sparsity is the sparsity bit pattern for z. On input it corresponds to \c G and on output it is undefined. */ template inline void reverse_sparse_jacobian_csum_op( size_t i_z , const addr_t* arg , Vector_set& sparsity ) { size_t i, j; i = arg[0] + arg[1]; j = 2; while(i--) { ++j; CPPAD_ASSERT_UNKNOWN( size_t(arg[j]) < i_z ); sparsity.binary_union( arg[j] , // index in sparsity for result arg[j] , // index in sparsity for left operand i_z , // index for right operand sparsity // sparsity vector for right operand ); } } /*! Reverse mode Hessian sparsity pattern for CSumOp operator. This operation is \verbatim z = q + x(1) + ... + x(m) - y(1) - ... - y(n). H(y, x, w, ...) = G[ z(x, y), y, x, w, ... ] \endverbatim \tparam Vector_set is the type used for vectors of sets. It can be either \c sparse_pack, \c sparse_set, or \c sparse_list. \param i_z variable index corresponding to the result for this operation; i.e. the index in \a sparsity corresponding to z. \param arg \a arg[0] is the number of addition variables in this cummulative summation; i.e., m + n. \n \a arg[1] is the number of subtraction variables in this cummulative summation; i.e., \c m. \n parameter[ arg[2] ] is the parameter value \c q in this cummunative summation. \n arg[2+i] for i = 1 , ... , m is the value x(i). \n arg[2+arg[0]+i] for i = 1 , ... , n is the value y(i). \param rev_jacobian rev_jacobian[i_z] is all false (true) if the Jabobian of G with respect to z must be zero (may be non-zero). \n \n For i = 1 , ... , m rev_jacobian[ arg[2+i] ] is all false (true) if the Jacobian with respect to x(i) is zero (may be non-zero). On input, it corresponds to the function G, and on output it corresponds to the function H. \n \n For i = 1 , ... , n rev_jacobian[ arg[2+arg[0]+i] ] is all false (true) if the Jacobian with respect to y(i) is zero (may be non-zero). On input, it corresponds to the function G, and on output it corresponds to the function H. \param rev_hes_sparsity The set with index \a i_z in in \a rev_hes_sparsity is the Hessian sparsity pattern for the fucntion G where one of the partials derivative is with respect to z. \n \n For i = 1 , ... , m The set with index arg[2+i] in \a rev_hes_sparsity is the Hessian sparsity pattern where one of the partials derivative is with respect to x(i). On input, it corresponds to the function G, and on output it corresponds to the function H. \n \n For i = 1 , ... , n The set with index arg[2+arg[0]+i] in \a rev_hes_sparsity is the Hessian sparsity pattern where one of the partials derivative is with respect to y(i). On input, it corresponds to the function G, and on output it corresponds to the function H. */ template inline void reverse_sparse_hessian_csum_op( size_t i_z , const addr_t* arg , bool* rev_jacobian , Vector_set& rev_hes_sparsity ) { size_t i, j; i = arg[0] + arg[1]; j = 2; while(i--) { ++j; CPPAD_ASSERT_UNKNOWN( size_t(arg[j]) < i_z ); rev_hes_sparsity.binary_union( arg[j] , // index in sparsity for result arg[j] , // index in sparsity for left operand i_z , // index for right operand rev_hes_sparsity // sparsity vector for right operand ); rev_jacobian[arg[j]] |= rev_jacobian[i_z]; } } } // END_CPPAD_NAMESPACE # endif