CppAD: A C++ Algorithmic Differentiation Package 20110419
for_jac_sweep.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_FOR_JAC_SWEEP_INCLUDED
00003 # define CPPAD_FOR_JAC_SWEEP_INCLUDED
00004 /* --------------------------------------------------------------------------
00005 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-11 Bradley M. Bell
00006 
00007 CppAD is distributed under multiple licenses. This distribution is under
00008 the terms of the 
00009                     Common Public License Version 1.0.
00010 
00011 A copy of this license is included in the COPYING file of this distribution.
00012 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00013 -------------------------------------------------------------------------- */
00014 # include <set>
00015 
00016 CPPAD_BEGIN_NAMESPACE
00017 /*!
00018 \file for_jac_sweep.hpp
00019 Compute Forward mode Jacobian sparsity patterns.
00020 */
00021 
00022 /*!
00023 \def CPPAD_FOR_JAC_SWEEP_TRACE
00024 This value is either zero or one. 
00025 Zero is the normal operational value.
00026 If it is one, a trace of every for_jac_sweep computation is printed.
00027 */
00028 # define CPPAD_FOR_JAC_SWEEP_TRACE 0
00029 
00030 /*!
00031 Given the sparsity pattern for the independent variables,
00032 ForJacSweep computes the sparsity pattern for all the other variables.
00033 
00034 \tparam Base
00035 base type for the operator; i.e., this operation sequence was recorded
00036 using AD< \a Base > and computations by this routine are done using type 
00037 \a Base.
00038 
00039 \tparam Vector_set
00040 is the type used for vectors of sets. It can be either
00041 \c sparse_pack or \c sparse_set.
00042 
00043 \param n
00044 is the number of independent variables on the tape.
00045 
00046 \param numvar
00047 is the total number of variables on the tape; i.e.,
00048 \a play->num_rec_var().
00049 
00050 \param play
00051 The information stored in \a play
00052 is a recording of the operations corresponding to a function
00053 \f[
00054         F : {\bf R}^n \rightarrow {\bf R}^m
00055 \f]
00056 where \f$ n \f$ is the number of independent variables
00057 and \f$ m \f$ is the number of dependent variables.
00058 The object \a play is effectly constant.
00059 It is not declared const because while playing back the tape
00060 the object \a play holds information about the currentl location
00061 with in the tape and this changes during playback.
00062 
00063 \param var_sparsity
00064 \b Input: For j = 1 , ... , \a n, 
00065 the sparsity pattern for the independent variable with index (j-1)
00066 corresponds to the set with index j in \a var_sparsity.
00067 \n
00068 \n
00069 \b Output: For i = \a n + 1 , ... , \a numvar - 1,
00070 the sparsity pattern for the variable with index i on the tape
00071 corresponds to the set with index i in \a var_sparsity.
00072 
00073 \par Checked Assertions:
00074 \li numvar == var_sparsity.n_set()
00075 \li numvar == play->num_rec_var()
00076 */
00077 
00078 template <class Base, class Vector_set>
00079 void ForJacSweep(
00080         size_t                n            ,
00081         size_t                numvar       ,
00082         player<Base>*         play         ,
00083         Vector_set&           var_sparsity )
00084 {
00085         OpCode           op;
00086         size_t         i_op;
00087         size_t        i_var;
00088 
00089         const size_t   *arg = 0;
00090         const size_t *arg_0 = 0;
00091 
00092         size_t            i, j, k;
00093 
00094         // check numvar argument
00095         CPPAD_ASSERT_UNKNOWN( play->num_rec_var()  == numvar );
00096         CPPAD_ASSERT_UNKNOWN( var_sparsity.n_set() == numvar );
00097 
00098         // length of the parameter vector (used by CppAD assert macros)
00099         const size_t num_par = play->num_rec_par();
00100 
00101         // cum_sparsity accumulates sparsity pattern a cummulative sum
00102         size_t limit = var_sparsity.end();
00103 
00104         // vecad_sparsity contains a sparsity pattern from each VecAD object
00105         // to all the other variables.
00106         // vecad_ind maps a VecAD index (the beginning of the
00107         // VecAD object) to its from index in vecad_sparsity
00108         size_t num_vecad_ind   = play->num_rec_vecad_ind();
00109         size_t num_vecad_vec   = play->num_rec_vecad_vec();
00110         Vector_set  vecad_sparsity;
00111         vecad_sparsity.resize(num_vecad_vec, limit);
00112         size_t* vecad_ind      = CPPAD_NULL;
00113         if( num_vecad_vec > 0 )
00114         {       size_t length;
00115                 vecad_ind     = CPPAD_TRACK_NEW_VEC(num_vecad_ind, vecad_ind);
00116                 j             = 0;
00117                 for(i = 0; i < num_vecad_vec; i++)
00118                 {       // length of this VecAD
00119                         length   = play->GetVecInd(j);
00120                         // set to proper index for this VecAD
00121                         vecad_ind[j] = i; 
00122                         for(k = 1; k <= length; k++)
00123                                 vecad_ind[j+k] = num_vecad_vec; // invalid index
00124                         // start of next VecAD
00125                         j       += length + 1;
00126                 }
00127                 CPPAD_ASSERT_UNKNOWN( j == play->num_rec_vecad_ind() );
00128         }
00129 
00130         // work space used by UserOp.
00131         typedef std::set<size_t> size_set;
00132         const size_t user_q = limit; // maximum element plus one
00133         size_set::iterator set_itr;  // iterator for a standard set
00134         size_set::iterator set_end;  // end of iterator sequence
00135         vector< size_set > user_r;   // sparsity pattern for the argument x
00136         vector< size_set > user_s;   // sparisty pattern for the result y
00137         size_t user_index = 0;       // indentifier for this user_atomic operation
00138         size_t user_id    = 0;       // user identifier for this call to operator
00139         size_t user_i     = 0;       // index in result vector
00140         size_t user_j     = 0;       // index in argument vector
00141         size_t user_m     = 0;       // size of result vector
00142         size_t user_n     = 0;       // size of arugment vector
00143         // next expected operator in a UserOp sequence
00144         enum { user_start, user_arg, user_ret, user_end } user_state = user_start;
00145 
00146 # if CPPAD_FOR_JAC_SWEEP_TRACE
00147         std::cout << std::endl;
00148         CppAD::vector<bool> z_value(limit);
00149 # endif
00150 
00151         // skip the BeginOp at the beginning of the recording
00152         play->start_forward(op, arg, i_op, i_var);
00153         CPPAD_ASSERT_UNKNOWN( op == BeginOp );
00154         arg_0 = arg;
00155         while(op != EndOp)
00156         {
00157                 // this op
00158                 play->next_forward(op, arg, i_op, i_var);
00159                 CPPAD_ASSERT_UNKNOWN( (i_op > n)  | (op == InvOp) );  
00160                 CPPAD_ASSERT_UNKNOWN( (i_op <= n) | (op != InvOp) );  
00161 
00162                 // rest of information depends on the case
00163                 switch( op )
00164                 {
00165                         case AbsOp:
00166                         CPPAD_ASSERT_NARG_NRES(op, 1, 1);
00167                         forward_sparse_jacobian_unary_op(
00168                                 i_var, arg[0], var_sparsity
00169                         );
00170                         break;
00171                         // -------------------------------------------------
00172 
00173                         case AddvvOp:
00174                         CPPAD_ASSERT_NARG_NRES(op, 2, 1);
00175                         forward_sparse_jacobian_binary_op(
00176                                 i_var, arg, var_sparsity
00177                         );
00178                         break;
00179                         // -------------------------------------------------
00180 
00181                         case AddpvOp:
00182                         CPPAD_ASSERT_NARG_NRES(op, 2, 1);
00183                         forward_sparse_jacobian_unary_op(
00184                                 i_var, arg[1], var_sparsity
00185                         );
00186                         break;
00187                         // -------------------------------------------------
00188 
00189                         case AcosOp:
00190                         // acos(x) and sqrt(1 - x * x) are computed in pairs
00191                         // but ivar + 1 should only be used here
00192                         CPPAD_ASSERT_NARG_NRES(op, 1, 2);
00193                         forward_sparse_jacobian_unary_op(
00194                                 i_var, arg[0], var_sparsity
00195                         );
00196                         break;
00197                         // -------------------------------------------------
00198 
00199                         case AsinOp:
00200                         // asin(x) and sqrt(1 - x * x) are computed in pairs
00201                         // but ivar + 1 should only be used here
00202                         CPPAD_ASSERT_NARG_NRES(op, 1, 2);
00203                         forward_sparse_jacobian_unary_op(
00204                                 i_var, arg[0], var_sparsity
00205                         );
00206                         break;
00207                         // -------------------------------------------------
00208 
00209                         case AtanOp:
00210                         // atan(x) and 1 + x * x must be computed in pairs
00211                         // but ivar + 1 should only be used here
00212                         CPPAD_ASSERT_NARG_NRES(op, 1, 2);
00213                         forward_sparse_jacobian_unary_op(
00214                                 i_var, arg[0], var_sparsity
00215                         );
00216                         break;
00217                         // -------------------------------------------------
00218 
00219                         case CSumOp:
00220                         // CSumOp has a variable number of arguments and
00221                         // next_forward thinks it one has one argument.
00222                         // we must inform next_forward of this special case.
00223                         play->forward_csum(op, arg, i_op, i_var);
00224                         forward_sparse_jacobian_csum_op(
00225                                 i_var, arg, var_sparsity
00226                         );
00227                         break;
00228                         // -------------------------------------------------
00229 
00230                         case CExpOp:
00231                         forward_sparse_jacobian_cond_op(
00232                                 i_var, arg, num_par, var_sparsity
00233                         );
00234                         break;
00235                         // ---------------------------------------------------
00236 
00237                         case ComOp:
00238                         CPPAD_ASSERT_NARG_NRES(op, 4, 0);
00239                         CPPAD_ASSERT_UNKNOWN( arg[1] > 1 );
00240                         break;
00241                         // --------------------------------------------------
00242 
00243                         case CosOp:
00244                         // cosine and sine must come in pairs
00245                         // but ivar should only be used here
00246                         CPPAD_ASSERT_NARG_NRES(op, 1, 2);
00247                         forward_sparse_jacobian_unary_op(
00248                                 i_var, arg[0], var_sparsity
00249                         );
00250                         break;
00251                         // ---------------------------------------------------
00252 
00253                         case CoshOp:
00254                         // hyperbolic cosine and sine must come in pairs
00255                         // but ivar should only be used here
00256                         CPPAD_ASSERT_NARG_NRES(op, 1, 2);
00257                         forward_sparse_jacobian_unary_op(
00258                                 i_var, arg[0], var_sparsity
00259                         );
00260                         break;
00261                         // -------------------------------------------------
00262 
00263                         case DisOp:
00264                         CPPAD_ASSERT_NARG_NRES(op, 2, 1);
00265                         var_sparsity.clear(i_var);
00266                         break;
00267                         // -------------------------------------------------
00268 
00269                         case DivvvOp:
00270                         CPPAD_ASSERT_NARG_NRES(op, 2, 1);
00271                         forward_sparse_jacobian_binary_op(
00272                                 i_var, arg, var_sparsity
00273                         );
00274                         break;
00275                         // -------------------------------------------------
00276 
00277                         case DivpvOp:
00278                         CPPAD_ASSERT_NARG_NRES(op, 2, 1);
00279                         forward_sparse_jacobian_unary_op(
00280                                 i_var, arg[1], var_sparsity
00281                         );
00282                         break;
00283                         // -------------------------------------------------
00284 
00285                         case DivvpOp:
00286                         CPPAD_ASSERT_NARG_NRES(op, 2, 1);
00287                         forward_sparse_jacobian_unary_op(
00288                                 i_var, arg[0], var_sparsity
00289                         );
00290                         break;
00291                         // -------------------------------------------------
00292 
00293                         case EndOp:
00294                         CPPAD_ASSERT_NARG_NRES(op, 0, 0);
00295                         break;
00296                         // -------------------------------------------------
00297 
00298                         case ExpOp:
00299                         CPPAD_ASSERT_NARG_NRES(op, 1, 1);
00300                         forward_sparse_jacobian_unary_op(
00301                                 i_var, arg[0], var_sparsity
00302                         );
00303                         break;
00304                         // -------------------------------------------------
00305 
00306                         case InvOp:
00307                         CPPAD_ASSERT_NARG_NRES(op, 0, 1);
00308                         // sparsity pattern is already defined
00309                         break;
00310                         // -------------------------------------------------
00311 
00312                         case LdpOp:
00313                         forward_sparse_load_op(
00314                                 op,
00315                                 i_var,
00316                                 arg,
00317                                 num_vecad_ind,
00318                                 vecad_ind,
00319                                 var_sparsity,
00320                                 vecad_sparsity
00321                         );
00322                         break;
00323                         // -------------------------------------------------
00324 
00325                         case LdvOp:
00326                         forward_sparse_load_op(
00327                                 op,
00328                                 i_var,
00329                                 arg,
00330                                 num_vecad_ind,
00331                                 vecad_ind,
00332                                 var_sparsity,
00333                                 vecad_sparsity
00334                         );
00335                         break;
00336                         // -------------------------------------------------
00337 
00338                         case LogOp:
00339                         CPPAD_ASSERT_NARG_NRES(op, 1, 1);
00340                         forward_sparse_jacobian_unary_op(
00341                                 i_var, arg[0], var_sparsity
00342                         );
00343                         break;
00344                         // -------------------------------------------------
00345 
00346                         case MulvvOp:
00347                         CPPAD_ASSERT_NARG_NRES(op, 2, 1);
00348                         forward_sparse_jacobian_binary_op(
00349                                 i_var, arg, var_sparsity
00350                         );
00351                         break;
00352                         // -------------------------------------------------
00353 
00354                         case MulpvOp:
00355                         CPPAD_ASSERT_NARG_NRES(op, 2, 1);
00356                         forward_sparse_jacobian_unary_op(
00357                                 i_var, arg[1], var_sparsity
00358                         );
00359                         break;
00360                         // -------------------------------------------------
00361 
00362                         case ParOp:
00363                         CPPAD_ASSERT_NARG_NRES(op, 1, 1);
00364                         var_sparsity.clear(i_var);
00365                         break;
00366                         // -------------------------------------------------
00367 
00368                         case PowvpOp:
00369                         CPPAD_ASSERT_NARG_NRES(op, 2, 3);
00370                         forward_sparse_jacobian_unary_op(
00371                                 i_var, arg[0], var_sparsity
00372                         );
00373                         break;
00374                         // -------------------------------------------------
00375 
00376                         case PowpvOp:
00377                         CPPAD_ASSERT_NARG_NRES(op, 2, 3);
00378                         forward_sparse_jacobian_unary_op(
00379                                 i_var, arg[1], var_sparsity
00380                         );
00381                         break;
00382                         // -------------------------------------------------
00383 
00384                         case PowvvOp:
00385                         CPPAD_ASSERT_NARG_NRES(op, 2, 3);
00386                         forward_sparse_jacobian_binary_op(
00387                                 i_var, arg, var_sparsity
00388                         );
00389                         break;
00390                         // -------------------------------------------------
00391 
00392                         case PripOp:
00393                         CPPAD_ASSERT_NARG_NRES(op, 2, 0);
00394                         var_sparsity.clear(i_var);
00395                         break;
00396                         // -------------------------------------------------
00397 
00398                         case PrivOp:
00399                         CPPAD_ASSERT_NARG_NRES(op, 2, 0);
00400                         break;
00401                         // -------------------------------------------------
00402 
00403                         case SinOp:
00404                         // sine and cosine must come in pairs
00405                         // but ivar should only be used here
00406                         CPPAD_ASSERT_NARG_NRES(op, 1, 2);
00407                         forward_sparse_jacobian_unary_op(
00408                                 i_var, arg[0], var_sparsity
00409                         );
00410                         break;
00411                         // -------------------------------------------------
00412 
00413                         case SinhOp:
00414                         // hyperbolic sine and cosine must come in pairs
00415                         // but ivar should only be used here
00416                         CPPAD_ASSERT_NARG_NRES(op, 1, 2);
00417                         forward_sparse_jacobian_unary_op(
00418                                 i_var, arg[0], var_sparsity
00419                         );
00420                         break;
00421                         // -------------------------------------------------
00422 
00423                         case SqrtOp:
00424                         CPPAD_ASSERT_NARG_NRES(op, 1, 1);
00425                         forward_sparse_jacobian_unary_op(
00426                                 i_var, arg[0], var_sparsity
00427                         );
00428                         break;
00429                         // -------------------------------------------------
00430 
00431                         case StppOp:
00432                         CPPAD_ASSERT_NARG_NRES(op, 3, 0);
00433                         // storing a parameter does not affect vector sparsity
00434                         break;
00435                         // -------------------------------------------------
00436 
00437                         case StpvOp:
00438                         forward_sparse_store_op(
00439                                 op,
00440                                 arg,
00441                                 num_vecad_ind,
00442                                 vecad_ind,
00443                                 var_sparsity,
00444                                 vecad_sparsity
00445                         );
00446                         break;
00447                         // -------------------------------------------------
00448 
00449                         case StvpOp:
00450                         CPPAD_ASSERT_NARG_NRES(op, 3, 0);
00451                         // storing a parameter does not affect vector sparsity
00452                         break;
00453                         // -------------------------------------------------
00454 
00455                         case StvvOp:
00456                         forward_sparse_store_op(
00457                                 op,
00458                                 arg,
00459                                 num_vecad_ind,
00460                                 vecad_ind,
00461                                 var_sparsity,
00462                                 vecad_sparsity
00463                         );
00464                         break;
00465                         // -------------------------------------------------
00466 
00467                         case SubvvOp:
00468                         CPPAD_ASSERT_NARG_NRES(op, 2, 1);
00469                         forward_sparse_jacobian_binary_op(
00470                                 i_var, arg, var_sparsity
00471                         );
00472                         break;
00473                         // -------------------------------------------------
00474 
00475                         case SubpvOp:
00476                         CPPAD_ASSERT_NARG_NRES(op, 2, 1);
00477                         forward_sparse_jacobian_unary_op(
00478                                 i_var, arg[1], var_sparsity
00479                         );
00480                         break;
00481                         // -------------------------------------------------
00482 
00483                         case SubvpOp:
00484                         CPPAD_ASSERT_NARG_NRES(op, 2, 1);
00485                         forward_sparse_jacobian_unary_op(
00486                                 i_var, arg[0], var_sparsity
00487                         );
00488                         break;
00489                         // -------------------------------------------------
00490 
00491                         case UserOp:
00492                         // start or end an atomic operation sequence
00493                         CPPAD_ASSERT_UNKNOWN( NumRes( UserOp ) == 0 );
00494                         CPPAD_ASSERT_UNKNOWN( NumArg( UserOp ) == 4 );
00495                         if( user_state == user_start )
00496                         {       user_index = arg[0];
00497                                 user_id    = arg[1];
00498                                 user_n     = arg[2];
00499                                 user_m     = arg[3];
00500                                 if(user_r.size() < user_n )
00501                                         user_r.resize(user_n);
00502                                 if(user_s.size() < user_m )
00503                                         user_s.resize(user_m);
00504                                 user_j     = 0;
00505                                 user_i     = 0;
00506                                 user_state = user_arg;
00507                         }
00508                         else
00509                         {       CPPAD_ASSERT_UNKNOWN( user_state == user_end );
00510                                 CPPAD_ASSERT_UNKNOWN( user_index == arg[0] );
00511                                 CPPAD_ASSERT_UNKNOWN( user_id    == arg[1] );
00512                                 CPPAD_ASSERT_UNKNOWN( user_n     == arg[2] );
00513                                 CPPAD_ASSERT_UNKNOWN( user_m     == arg[3] );
00514                                 user_state = user_start;
00515                         }
00516                         break;
00517 
00518                         case UsrapOp:
00519                         // parameter argument in an atomic operation sequence
00520                         CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
00521                         CPPAD_ASSERT_UNKNOWN( user_j < user_n );
00522                         CPPAD_ASSERT_UNKNOWN( arg[0] < num_par );
00523                         // parameters have an empty sparsity pattern
00524                         user_r[user_j].clear();
00525                         ++user_j;
00526                         if( user_j == user_n )
00527                         {       // call users function for this operation
00528                                 user_atomic<Base>::for_jac_sparse(user_index, user_id,
00529                                         user_n, user_m, user_q, user_r, user_s
00530                                 );
00531                                 user_state = user_ret;
00532                         }
00533                         break;
00534 
00535                         case UsravOp:
00536                         // variable argument in an atomic operation sequence
00537                         CPPAD_ASSERT_UNKNOWN( user_state == user_arg );
00538                         CPPAD_ASSERT_UNKNOWN( user_j < user_n );
00539                         CPPAD_ASSERT_UNKNOWN( arg[0] <= i_var );
00540                         // set user_r[user_j] to sparsity pattern for variable arg[0]
00541                         user_r[user_j].clear();
00542                         var_sparsity.begin(arg[0]);
00543                         i = var_sparsity.next_element();
00544                         while( i < user_q )
00545                         {       user_r[user_j].insert(i);
00546                                 i = var_sparsity.next_element();
00547                         }
00548                         ++user_j;
00549                         if( user_j == user_n )
00550                         {       // call users function for this operation
00551                                 user_atomic<Base>::for_jac_sparse(user_index, user_id,
00552                                         user_n, user_m, user_q, user_r, user_s
00553                                 );
00554                                 user_state = user_ret;
00555                         }
00556                         break;
00557 
00558                         case UsrrpOp:
00559                         // parameter result in an atomic operation sequence
00560                         CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
00561                         CPPAD_ASSERT_UNKNOWN( user_i < user_m );
00562                         user_i++;
00563                         if( user_i == user_m )
00564                                 user_state = user_end;
00565                         break;
00566 
00567                         case UsrrvOp:
00568                         // variable result in an atomic operation sequence
00569                         CPPAD_ASSERT_UNKNOWN( user_state == user_ret );
00570                         CPPAD_ASSERT_UNKNOWN( user_i < user_m );
00571                         // It might be faster if we add set union to var_sparsity
00572                         // where one of the sets is not in var_sparsity
00573                         set_itr = user_s[user_i].begin();
00574                         set_end = user_s[user_i].end();
00575                         while( set_itr != set_end )
00576                                 var_sparsity.add_element(i_var, *set_itr++);
00577                         user_i++;
00578                         if( user_i == user_m )
00579                                 user_state = user_end;
00580                         break;
00581                         // -------------------------------------------------
00582 
00583                         default:
00584                         CPPAD_ASSERT_UNKNOWN(0);
00585                 }
00586 # if CPPAD_FOR_JAC_SWEEP_TRACE
00587                 // value for this variable
00588                 for(j = 0; j < limit; j++)
00589                         z_value[j] = false;
00590                 var_sparsity.begin(i_var);
00591                 j = var_sparsity.next_element();
00592                 while( j < limit )
00593                 {       z_value[j] = true;
00594                         j = var_sparsity.next_element();
00595                 }
00596                 printOp(
00597                         std::cout,
00598                         play,
00599                         i_var,
00600                         op,
00601                         arg,
00602                         1,
00603                         &z_value,
00604                         0,
00605                         (CppAD::vector<bool> *) CPPAD_NULL
00606                 );
00607         }
00608         std::cout << std::endl;
00609 # else
00610         }
00611 # endif
00612         CPPAD_ASSERT_UNKNOWN( i_var + 1 == play->num_rec_var() );
00613 
00614         if( vecad_ind != CPPAD_NULL )
00615                 CPPAD_TRACK_DEL_VEC(vecad_ind);
00616 
00617         return;
00618 }
00619 
00620 # undef CPPAD_FOR_JAC_SWEEP_TRACE
00621 
00622 CPPAD_END_NAMESPACE
00623 # endif