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