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