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