Changeset 1469


Ignore:
Timestamp:
Jul 30, 2009 10:45:53 AM (11 years ago)
Author:
bradbell
Message:

trunk: Add reverse sparse Jacobian and Hessian CondExp? routines.

test_more.cpp: change RevSparseJac? to rev_sparse_jac.
cond_exp_ad.cpp: new pattern valid for all independent variable values.
rev_sparse_jac.cpp: new pattern valid for all independent variable values.
rev_sparse_jac.cpp: new pattern valid for all independent variable values.
for_jac_sweep.hpp: add some white space.
rev_jac_sweep.hpp: use new sparse Jacobian routine.
rev_hes_sweep.hpp: use new sparse Hessian routine.
cond_op.hpp: define new routines.

Location:
trunk
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • trunk/cppad/local/cond_op.hpp

    r1468 r1469  
    407407}
    408408
     409/*!
     410Compute reverse Jacobian sparsity patterns for op = CExpOp.
     411
     412This routine is given the sparsity patterns
     413for a function G(z, y, x, ... )
     414and it uses them to compute the sparsity patterns for
     415\verbatim
     416        H( y, x, w , u , ... ) = G[ z(x,y) , y , x , w , u , ... ]
     417\endverbatim
     418where y represents the combination of y_0, y_1, y_2, and y_3.
     419
     420\copydetails sparse_conditional_exp_op
     421
     422\param sparsity
     423If y_2 is a variable,
     424for k = 0 , ... , nc_sparsity-1,
     425\a sparsity [ \a arg[4] * nc_sparsity + k ]
     426is the sparsity bit pattern corresponding to y_2.
     427On input, this pattern corresponds to the function G.
     428On ouput, it corresponds to the function H.
     429\n
     430\n
     431If y_3 is a variable,
     432for k = 0 , ... , nc_sparsity-1,
     433\a sparsity [ \a arg[5] * nc_sparsity + k ]
     434is the sparsity bit pattern corresponding to y_3.
     435On input, this pattern corresponds to the function G.
     436On ouput, it corresponds to the function H.
     437\n
     438\n
     439For k = 0 , ... , nc_sparsity-1,
     440\a sparsity [ \a i_z * nc_sparsity + k ]
     441is the sparsity bit pattern corresponding to z.
     442On input and output, this pattern corresponds to the function G.
     443*/
     444template <class Pack>
     445inline void reverse_sparse_jacobian_cond_op(
     446        size_t         i_z           ,
     447        const size_t*  arg           ,
     448        size_t         num_par       ,
     449        size_t         nc_sparsity   ,
     450        Pack*          sparsity      )
     451{       
     452        Pack* z   = sparsity + i_z * nc_sparsity;
     453
     454        CPPAD_ASSERT_UNKNOWN( arg[0] < static_cast<size_t> (CompareNe) );
     455        CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 );
     456        CPPAD_ASSERT_UNKNOWN( NumRes(CExpOp) == 1 );
     457        CPPAD_ASSERT_UNKNOWN( arg[1] != 0 );
     458
     459# ifndef NDEBUG
     460        if( arg[1] & 1 )
     461        {       CPPAD_ASSERT_UNKNOWN( arg[2] < i_z );
     462        }
     463        else
     464        {       CPPAD_ASSERT_UNKNOWN( arg[2] < num_par );
     465        }
     466        if( arg[1] & 2 )
     467        {       CPPAD_ASSERT_UNKNOWN( arg[3] < i_z );
     468        }
     469        else
     470        {       CPPAD_ASSERT_UNKNOWN( arg[3] < num_par );
     471        }
     472        if( ! ( arg[1] & 4 ) )
     473        {       CPPAD_ASSERT_UNKNOWN( arg[4] < num_par );
     474        }
     475        if( ! ( arg[1] & 8 ) )
     476        {       CPPAD_ASSERT_UNKNOWN( arg[5] < num_par );
     477        }
     478# endif
     479        size_t k;
     480        if( arg[1] & 4 )
     481        {       CPPAD_ASSERT_UNKNOWN( arg[4] < i_z );
     482                Pack* y_2 = sparsity + arg[4] * nc_sparsity;
     483                k = nc_sparsity;
     484                while(k--)
     485                        y_2[k] |= z[k];
     486        }
     487        if( arg[1] & 8 )
     488        {       CPPAD_ASSERT_UNKNOWN( arg[5] < i_z );
     489                Pack* y_3 = sparsity + arg[5] * nc_sparsity;
     490                k = nc_sparsity;
     491                while(k--)
     492                        y_3[k] |= z[k];
     493        }
     494        return;
     495}
     496
     497/*!
     498Compute reverse Hessian sparsity patterns for op = CExpOp.
     499
     500This routine is given the sparsity patterns
     501for a function G(z, y, x, ... )
     502and it uses them to compute the sparsity patterns for
     503\verbatim
     504        H( y, x, w , u , ... ) = G[ z(x,y) , y , x , w , u , ... ]
     505\endverbatim
     506where y represents the combination of y_0, y_1, y_2, and y_3.
     507
     508\copydetails sparse_conditional_exp_op
     509
     510\param z_jac
     511is all true (ones complement of 0) if the scalar valued
     512function we are computing the Hessian sparsity for
     513has a non-zero partial with respect to the variable z
     514(actually may have a non-zero partial with respect to z).
     515Otherwise it zero.
     516
     517\param hes_sparsity
     518For k = 0 , ... , nc_sparsity-1,
     519\a sparsity [ \a i_z * nc_sparsity + k ]
     520is the sparsity bit pattern corresponding to z.
     521On input and output, this pattern corresponds to the function G.
     522\n
     523\n
     524For k = 0 , ... , nc_sparsity-1,
     525\a sparsity [ \a arg[4] * nc_sparsity + k ]
     526is the sparsity bit pattern corresponding to y_2.
     527On input, this pattern corresponds to the function G.
     528On output, this pattern corresponds to the function H.
     529\n
     530\n
     531For k = 0 , ... , nc_sparsity-1,
     532\a sparsity [ \a arg[5] * nc_sparsity + k ]
     533is the sparsity bit pattern corresponding to y_3.
     534On input, this pattern corresponds to the function G.
     535On output, this pattern corresponds to the function H.
     536*/
     537template <class Pack>
     538inline void reverse_sparse_hessian_cond_op(
     539        size_t         i_z           ,
     540        const size_t*  arg           ,
     541        size_t         num_par       ,
     542        Pack           z_jac         ,
     543        size_t         nc_sparsity   ,
     544        Pack*          hes_sparsity  )
     545{       
     546
     547        CPPAD_ASSERT_UNKNOWN( arg[0] < static_cast<size_t> (CompareNe) );
     548        CPPAD_ASSERT_UNKNOWN( NumArg(CExpOp) == 6 );
     549        CPPAD_ASSERT_UNKNOWN( NumRes(CExpOp) == 1 );
     550        CPPAD_ASSERT_UNKNOWN( arg[1] != 0 );
     551
     552        Pack* z_hes   = hes_sparsity + i_z * nc_sparsity;
     553
     554# ifndef NDEBUG
     555        if( arg[1] & 1 )
     556        {       CPPAD_ASSERT_UNKNOWN( arg[2] < i_z );
     557        }
     558        else
     559        {       CPPAD_ASSERT_UNKNOWN( arg[2] < num_par );
     560        }
     561        if( arg[1] & 2 )
     562        {       CPPAD_ASSERT_UNKNOWN( arg[3] < i_z );
     563        }
     564        else
     565        {       CPPAD_ASSERT_UNKNOWN( arg[3] < num_par );
     566        }
     567        if( ! ( arg[1] & 4 ) )
     568        {       CPPAD_ASSERT_UNKNOWN( arg[4] < num_par );
     569        }
     570        if( ! ( arg[1] & 8 ) )
     571        {       CPPAD_ASSERT_UNKNOWN( arg[5] < num_par );
     572        }
     573# endif
     574        size_t k;
     575        if( arg[1] & 4 )
     576        {       CPPAD_ASSERT_UNKNOWN( arg[4] < i_z );
     577                Pack* y_2_hes = hes_sparsity + arg[4] * nc_sparsity;
     578                k = nc_sparsity;
     579                while(k--)
     580                        y_2_hes[k] |= z_hes[k];
     581        }
     582        if( arg[1] & 8 )
     583        {       CPPAD_ASSERT_UNKNOWN( arg[5] < i_z );
     584                Pack* y_3_hes = hes_sparsity + arg[5] * nc_sparsity;
     585                k = nc_sparsity;
     586                while(k--)
     587                        y_3_hes[k] |= z_hes[k];
     588        }
     589        return;
     590}
     591
    409592CPPAD_END_NAMESPACE
    410593# endif
  • trunk/cppad/local/for_jac_sweep.hpp

    r1468 r1469  
    164164        // set the number of operators
    165165        const size_t numop_m1 = Rec->num_rec_op() - 1;
     166
    166167        // length of the parameter vector (used by CppAD assert macros)
    167168        const size_t num_par = Rec->num_rec_par();
  • trunk/cppad/local/rev_hes_sweep.hpp

    r1466 r1469  
    147147
    148148
    149         // used by CExp operator
    150         const Base *left = 0, *right = 0;
    151         Pack  *trueCaseh = 0;
    152         Pack  *falseCaseh = 0;
    153         Pack  zero(0);
     149        // length of the parameter vector (used by CppAD assert macros)
     150        const size_t num_par = Rec->num_rec_par();
    154151
    155152        size_t             j;
     
    262259                        // -------------------------------------------------
    263260                        case CExpOp:
    264                         CPPAD_ASSERT_UNKNOWN( n_res == 1);
    265                         CPPAD_ASSERT_UNKNOWN( n_arg == 6);
    266                         CPPAD_ASSERT_UNKNOWN( arg[1] != 0 );
    267 
    268                         if( arg[1] & 1 )
    269                                 left = Taylor + arg[2] * TaylorColDim;
    270                         else    left = Rec->GetPar(arg[2]);
    271                         if( arg[1] & 2 )
    272                                 right = Taylor + arg[3] * TaylorColDim;
    273                         else    right = Rec->GetPar(arg[3]);
    274                         if( arg[1] & 4 )
    275                         {       trueCaseh = RevHes + arg[4] * npv;
    276                                 for(j = 0; j < npv; j++)
    277                                 {       trueCaseh[j] |= CondExpTemplate(
    278                                                 CompareOp( arg[0] ),
    279                                                 *left,
    280                                                 *right,
    281                                                 Zh[j],
    282                                                 zero
    283                                         );
    284                                 }
    285                         }
    286                         if( arg[1] & 8 )
    287                         {       falseCaseh = RevHes + arg[5] * npv;
    288                                 for(j = 0; j < npv; j++)
    289                                 {       falseCaseh[j] |= CondExpTemplate(
    290                                                 CompareOp( arg[0] ),
    291                                                 *left,
    292                                                 *right,
    293                                                 zero,
    294                                                 Zh[j]
    295                                         );
    296                                 }
    297                         }
     261                        reverse_sparse_hessian_cond_op(
     262                                i_var, arg, num_par, *Zr, npv, RevHes
     263                        );
    298264                        break;
    299265                        // ---------------------------------------------------
  • trunk/cppad/local/rev_jac_sweep.hpp

    r1465 r1469  
    137137        size_t            j;
    138138
    139         // used by CExp operator
    140         const Base  *left = 0, *right = 0;
    141         Pack        *trueCase = 0, *falseCase = 0;
    142         Pack         zero(0);
     139        // length of the parameter vector (used by CppAD assert macros)
     140        const size_t num_par = Rec->num_rec_par();
    143141
    144142        // check numvar argument
     
    237235
    238236                        case CExpOp:
    239                         CPPAD_ASSERT_NARG_NRES(op, 6, 1);
    240                         CPPAD_ASSERT_UNKNOWN( arg[1] != 0 );
    241 
    242                         if( arg[1] & 1 )
    243                                 left = Taylor + arg[2] * TaylorColDim;
    244                         else    left = Rec->GetPar(arg[2]);
    245                         if( arg[1] & 2 )
    246                                 right = Taylor + arg[3] * TaylorColDim;
    247                         else    right = Rec->GetPar(arg[3]);
    248                         if( arg[1] & 4 )
    249                         {       trueCase = RevJac + arg[4] * npv;
    250                                 for(j = 0; j < npv; j++)
    251                                 {       trueCase[j] |= CondExpTemplate(
    252                                                 CompareOp( arg[0] ),
    253                                                 *left,
    254                                                 *right,
    255                                                 Z[j],
    256                                                 zero
    257                                         );
    258                                 }
    259                         }
    260                         if( arg[1] & 8 )
    261                         {       falseCase = RevJac + arg[5] * npv;
    262                                 for(j = 0; j < npv; j++)
    263                                 {       falseCase[j] |= CondExpTemplate(
    264                                                 CompareOp( arg[0] ),
    265                                                 *left,
    266                                                 *right,
    267                                                 zero,
    268                                                 Z[j]
    269                                         );
    270                                 }
    271                         }
     237                        reverse_sparse_jacobian_cond_op(
     238                                i_var, arg, num_par, npv, RevJac
     239                        );
    272240                        break;
    273241                        // ---------------------------------------------------
  • trunk/test_more/cond_exp_ad.cpp

    r1468 r1469  
    162162        for(i = 0; i < m; i++)
    163163        {       for(j = 0; j < n; j++)
    164                         ok &= ( Px[i * n + j] == (J[i * n + j] == 1.) );
     164                        ok &= ( Px[i * n + j] == ( j > 0 ) );
    165165        }
    166166
     
    315315        for(i = 0; i < m; i++)
    316316        {       for(j = 0; j < n; j++)
    317                         ok &= ( Px[i * n + j] == (J[i * n + j] == 1.) );
     317                        ok &= ( Px[i * n + j] == (j > 0) );
    318318        }
    319319
  • trunk/test_more/rev_sparse_hes.cpp

    r1370 r1469  
    11/* $Id$ */
    22/* --------------------------------------------------------------------------
    3 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell
     3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 Bradley M. Bell
    44
    55CppAD is distributed under multiple licenses. This distribution is under
     
    1313
    1414# include <cppad/cppad.hpp>
     15 
     16namespace { // Begin empty namespace
    1517
    16 
    17 bool RevSparseHes(void)
     18bool case_one(void)
    1819{       bool ok = true;
    1920        using namespace CppAD;
     
    5253        // CondExpLt(variable, variable, variable, variable)
    5354        sum += CondExpLt(X[1], X[2], sin(X[6]), cos(X[7]) );
    54         if( X[1] < X[2] )
    55                 Check[6 * n + 6] = true;
    56         else    Check[7 * n + 7] = true;
     55        Check[6 * n + 6] = true;
     56        Check[7 * n + 7] = true;
    5757       
    5858        // pow(variable, variable)
     
    102102        return ok;
    103103}
     104
     105bool case_two(void)
     106{       bool ok = true;
     107        using namespace CppAD;
     108
     109        // dimension of the domain space
     110        size_t n = 4;
     111
     112        // dimension of the range space
     113        size_t m = 1;
     114
     115        // temporary indices
     116        size_t i, j;
     117
     118        // initialize check values to false
     119        CPPAD_TEST_VECTOR<bool> Check(n * n);
     120        for(j = 0; j < n * n; j++)
     121                Check[j] = false;
     122
     123        // independent variable vector
     124        CPPAD_TEST_VECTOR< AD<double> > X(n);
     125        for(j = 0; j < n; j++)
     126                X[j] = AD<double>(j);
     127        Independent(X);
     128
     129        // Test the case where dependent variable is a non-linear function
     130        // of the result of a conditional expression.
     131        CPPAD_TEST_VECTOR< AD<double> > Y(m);
     132        Y[0] = CondExpLt(X[0], X[1], X[2], X[3]);
     133        Y[0] = cos(Y[0]) + X[0] + X[1];
     134
     135        // Hessian with respect to x[0] and x[1] is zero.
     136        // Hessian with respect to x[2] and x[3] is full
     137        // (although we know that there are no cross terms, this is an
     138        // inefficiency of the conditional expression operator).
     139        Check[2 * n + 2] = Check[ 2 * n + 3 ] = true;
     140        Check[3 * n + 2] = Check[ 3 * n + 3 ] = true;
     141       
     142        // create function object F : X -> Y
     143        ADFun<double> F(X, Y);
     144
     145        // sparsity pattern for the identity function U(x) = x
     146        CPPAD_TEST_VECTOR<bool> Px(n * n);
     147        for(i = 0; i < n; i++)
     148        {       for(j = 0; j < n; j++)
     149                        Px[ i * n + j ] = false;
     150                Px[ i * n + i ] = true;
     151        }
     152
     153        // compute sparsity pattern for Jacobian of F(U(x))
     154        F.ForSparseJac(n, Px);
     155
     156        // compute sparsity pattern for Hessian of F_0 ( U(x) )
     157        CPPAD_TEST_VECTOR<bool> Py(m);
     158        Py[0] = true;
     159        CPPAD_TEST_VECTOR<bool> Pxx(n * n);
     160        Pxx = F.RevSparseHes(n, Py);
     161
     162        // check values
     163        for(j = 0; j < n * n; j++)
     164                ok &= (Pxx[j] == Check[j]);
     165
     166        return ok;
     167}
     168
     169} // End of empty namespace
     170
     171bool rev_sparse_hes(void)
     172{       bool ok = true;
     173        // ok &= case_one();
     174        ok &= case_two();
     175
     176        return ok;
     177}
  • trunk/test_more/rev_sparse_jac.cpp

    r1370 r1469  
    11/* $Id$ */
    22/* --------------------------------------------------------------------------
    3 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell
     3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 Bradley M. Bell
    44
    55CppAD is distributed under multiple licenses. This distribution is under
     
    125125        Y[index] = CondExpLt(X[0], X[1], X[0], X[1]);
    126126        Check[index * n + 0] = true;
    127         Check[index * n + 1] = false;
     127        Check[index * n + 1] = true;
    128128        Check[index * n + 2] = false;
    129129        index++;
    130130        Y[index] = CondExpLt(X[0], X[1], AD<double>(3.), X[1]);
    131131        Check[index * n + 0] = false;
    132         Check[index * n + 1] = false;
     132        Check[index * n + 1] = true;
    133133        Check[index * n + 2] = false;
    134134        index++;
  • trunk/test_more/test_more.cpp

    r1401 r1469  
    6060extern bool PowInt(void);
    6161extern bool Reverse(void);
    62 extern bool RevSparseHes(void);
     62extern bool rev_sparse_hes(void);
    6363extern bool RevSparseJac(void);
    6464extern bool RevTwo(void);
     
    153153        ok &= Run( PowInt,          "PowInt"         );
    154154        ok &= Run( Reverse,         "Reverse"        );
    155         ok &= Run( RevSparseHes,    "RevSparseHes"  );
     155        ok &= Run( rev_sparse_hes,  "rev_sparse_hes" );
    156156        ok &= Run( RevSparseJac,    "RevSparseJac"   );
    157157        ok &= Run( RevTwo,          "RevTwo"         );
Note: See TracChangeset for help on using the changeset viewer.