Changeset 3712


Ignore:
Timestamp:
Aug 26, 2015 12:04:34 PM (5 years ago)
Author:
bradbell
Message:

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: 9068654d6c9e0157639951dd93dad71eb8149e77
end hash code: fc86286a32cd87ded0475c12f6a8ea5388bc0ccd

commit fc86286a32cd87ded0475c12f6a8ea5388bc0ccd
Author: Brad Bell <bradbell@…>
Date: Wed Aug 26 08:52:02 2015 -0700

Fix a bug in checkpoint revese mode Jacobian sparsity patterns for a subset of rows.
checkpoint.hpp: fix bug.
checkpoint.cpp: test for bug.

commit b7d27ce9be62bf303150df355fc369e1e668ff3f
Author: Brad Bell <bradbell@…>
Date: Tue Aug 25 17:55:36 2015 -0700

Use CPPAD_INTERNAL_SPARES_SET type for checkpoint sparsity information.

commit 9fa2f6bd5ac4a6541ca6d06b53ded266a46fcebb
Author: Brad Bell <bradbell@…>
Date: Sat Aug 22 18:09:13 2015 -0700

checkpoint.hpp: put calculation of sparsity in private function.
rev_sparse_jac.hpp: comments to separate routines.

commit 9d7fdd3ca29a7ccf2f02b25c07788af133284c04
Author: Brad Bell <bradbell@…>
Date: Sat Aug 22 17:05:02 2015 -0700

Change checkpoint sparse hessian to use internal pattern (less memory).

Location:
trunk
Files:
7 edited

Legend:

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

    r3709 r3712  
    326326                size_t q, const VectorSet &s, bool transpose = false
    327327        );
     328        // internal set sparsity version of RevSparseHes
     329        // (used by checkpoint functions only)
     330        void RevSparseHesCheckpoint(
     331                size_t                        q         ,
     332                vector<bool>&                 s         ,
     333                bool                          transpose ,
     334                CPPAD_INTERNAL_SPARSE_SET&    h
     335        );
     336        // internal set sparsity version of RevSparseJac
     337        // (used by checkpoint functions only)
     338        void RevSparseJacCheckpoint(
     339                size_t                        q          ,
     340                CPPAD_INTERNAL_SPARSE_SET&    r          ,
     341                bool                          transpose  ,
     342                bool                          dependency ,
     343                CPPAD_INTERNAL_SPARSE_SET&    s
     344        );
     345    // internal set sparsity version of RevSparseJac
     346    // (used by checkpoint functions only)
     347        void ForSparseJacCheckpoint(
     348        size_t                        q          ,
     349        CPPAD_INTERNAL_SPARSE_SET&    r          ,
     350        bool                          transpose  ,
     351        bool                          dependency ,
     352        CPPAD_INTERNAL_SPARSE_SET&    s
     353        );
    328354
    329355        /// amount of memory used for Jacobain sparsity pattern
  • trunk/cppad/local/checkpoint.hpp

    r3710 r3712  
    1313Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
    1414-------------------------------------------------------------------------- */
     15# include <cppad/local/sparse_set.hpp>
     16# include <cppad/local/sparse_list.hpp>
    1517
    1618namespace CppAD { // BEGIN_CPPAD_NAMESPACE
     
    178180template <class Base>
    179181class checkpoint : public atomic_base<Base> {
     182// ---------------------------------------------------------------------------
    180183private:
    181184        ADFun<Base> f_;
    182185        //
    183         // sparsity for f(x)^{(1)} (set by constructor)
    184         vector< std::set<size_t> > entire_jac_sparse_;
     186        /// sparsity for f(x)^{(1)}
     187        CPPAD_INTERNAL_SPARSE_SET entire_jac_sparse_;
    185188        //
    186         // sparsity for sum_i f_i(x)^{(2)}
    187         vector< std::set<size_t> > entire_hes_sparse_;
     189        /// sparsity for sum_i f_i(x)^{(2)}
     190        CPPAD_INTERNAL_SPARSE_SET  entire_hes_sparse_;
     191        //
     192        /// set entire_jac_sparse_
     193        void set_entire_jac_sparse(void)
     194        {       assert( entire_jac_sparse_.n_set() == 0 );
     195                bool transpose  = false;
     196                bool dependency = true;
     197                size_t n = f_.Domain();
     198                size_t m = f_.Range();
     199                // Use the choice for forward / reverse that results in smaller
     200                // size for the sparsity pattern of all variables in the tape.
     201                if( n <= m )
     202                {       CPPAD_INTERNAL_SPARSE_SET identity;
     203                        identity.resize(n, n);
     204                        for(size_t j = 0; j < n; j++)
     205                                identity.add_element(j, j);
     206                        f_.ForSparseJacCheckpoint(
     207                                n, identity, transpose, dependency, entire_jac_sparse_
     208                        );
     209                        f_.size_forward_set(0);
     210                }
     211                else
     212                {       CPPAD_INTERNAL_SPARSE_SET identity;
     213                        identity.resize(m, m);
     214                        for(size_t i = 0; i < m; i++)
     215                                identity.add_element(i, i);
     216                        f_.RevSparseJacCheckpoint(
     217                                m, identity, transpose, dependency, entire_jac_sparse_
     218                        );
     219                }
     220                CPPAD_ASSERT_UNKNOWN( f_.size_forward_set() == 0 );
     221                CPPAD_ASSERT_UNKNOWN( f_.size_forward_bool() == 0 );
     222        }
     223        //
     224        /// set entire_hes_sparse_
     225        void set_entire_hes_sparse(void)
     226        {       assert( entire_hes_sparse_.n_set() == 0 );
     227                size_t n = f_.Domain();
     228                size_t m = f_.Range();
     229                //
     230                // set version of sparsity for vector of all ones
     231                vector<bool> all_one(m);
     232                for(size_t i = 0; i < m; i++)
     233                        all_one[i] = true;
     234
     235                // set version of sparsity for n by n idendity matrix
     236                vector< std::set<size_t> > identity(n);
     237                for(size_t j = 0; j < n; j++)
     238                        identity[j].insert(j);
     239
     240                // compute sparsity pattern for H(x) = sum_i f_i(x)^{(2)}
     241                bool transpose  = false;
     242                bool dependency = false;
     243                f_.ForSparseJac(n, identity, transpose, dependency);
     244                f_.RevSparseHesCheckpoint(
     245                        n, all_one, transpose, entire_hes_sparse_
     246                );
     247                CPPAD_ASSERT_UNKNOWN( entire_hes_sparse_.n_set() == n );
     248                CPPAD_ASSERT_UNKNOWN( entire_hes_sparse_.end()   == n );
     249                //
     250                // drop the forward sparsity results from f_
     251                f_.size_forward_set(0);
     252        }
    188253public:
    189254        /*!
     
    227292                //
    228293                // set sparsity for entire Jacobian once and for all
    229                 assert( entire_jac_sparse_.size() == 0 );
    230                 bool transpose  = false;
    231                 bool dependency = true;
    232                 size_t n = f_.Domain();
    233                 size_t m = f_.Range();
    234                 // It is not clear if forward or reverse is best in sparse case,
    235                 // so use the best choice for the dense case (which this may be).
    236                 if( n <= m )
    237                 {       vector< std::set<size_t> > identity(n);
    238                         for(size_t j = 0; j < n; j++)
    239                                 identity[j].insert(j);
    240                         entire_jac_sparse_ = f_.ForSparseJac(
    241                                 n, identity, transpose, dependency
    242                         );
    243                         // drop the forward sparsity results from f_
    244                         f_.size_forward_set(0);
    245                 }
    246                 else
    247                 {       vector< std::set<size_t> > identity(m);
    248                         for(size_t i = 0; i < m; i++)
    249                                 identity[i].insert(i);
    250                         entire_jac_sparse_ = f_.RevSparseJac(
    251                                 m, identity, transpose, dependency
    252                         );
    253                 }
    254                 CPPAD_ASSERT_UNKNOWN( f_.size_forward_set() == 0 );
    255                 CPPAD_ASSERT_UNKNOWN( f_.size_forward_bool() == 0 );
     294                set_entire_jac_sparse();
    256295        }
    257296        /*!
     
    297336                const vector<Base>&      tx ,
    298337                      vector<Base>&      ty )
    299         {
     338        {       size_t n = f_.Domain();
     339                size_t m = f_.Range();
     340                //
    300341                CPPAD_ASSERT_UNKNOWN( f_.size_var() > 0 );
    301342                CPPAD_ASSERT_UNKNOWN( tx.size() % (q+1) == 0 );
    302343                CPPAD_ASSERT_UNKNOWN( ty.size() % (q+1) == 0 );
    303 # ifndef NDEBUG
    304                 size_t n = tx.size() / (q+1);
    305 # endif
    306                 size_t m = ty.size() / (q+1);
     344                CPPAD_ASSERT_UNKNOWN( n == tx.size() / (q+1) );
     345                CPPAD_ASSERT_UNKNOWN( m == ty.size() / (q+1) );
    307346                bool ok  = true;
    308                 size_t i, j;
    309 
    310                 // 2DO: test both forward and reverse vy information
     347                //
     348                // repeat it every time forward zero is used.
    311349                if( vx.size() > 0 )
    312350                {       // Compute Jacobian sparsity pattern.
    313                         assert( entire_jac_sparse_.size() > 0 );
    314                         for(i = 0; i < m; i++)
     351                        assert( entire_jac_sparse_.n_set() == m );
     352                        assert( entire_jac_sparse_.end()   == n );
     353                        for(size_t i = 0; i < m; i++)
    315354                        {       vy[i] = false;
    316                                 std::set<size_t>::const_iterator itr;
    317                                 const std::set<size_t>& s_i( entire_jac_sparse_[i] );
    318                                 for(itr = s_i.begin(); itr != s_i.end(); itr++)
    319                                 {       j = *itr;
    320                                         assert( j < n );
    321                                         // y[i] depends on the value of x[j]
     355                                entire_jac_sparse_.begin(i);
     356                                size_t j = entire_jac_sparse_.next_element();
     357                                while(j < n )
     358                                {       // y[i] depends on the value of x[j]
    322359                                        vy[i] |= vx[j];
     360                                        j = entire_jac_sparse_.next_element();
    323361                                }
    324362                        }
     
    388426                const vector< std::set<size_t> >&       r  ,
    389427                      vector< std::set<size_t> >&       s  )
    390         {       assert( entire_jac_sparse_.size() != 0 );
     428        {       assert( entire_jac_sparse_.n_set() != 0 );
    391429                assert( r.size() == f_.Domain() );
    392430                assert( s.size() == f_.Range() );
     
    394432                bool ok = true;
    395433                size_t m = f_.Range();
     434                size_t n = f_.Domain();
    396435                for(size_t i = 0; i < m; i++)
    397436                        s[i].clear();
     
    400439                for(size_t i = 0; i < m; i++)
    401440                {       // compute row i of the return pattern
    402                         std::set<size_t>::const_iterator itr_i;
    403                         const std::set<size_t>& jac_i( entire_jac_sparse_[i] );
    404                         for(itr_i = jac_i.begin(); itr_i != jac_i.end(); itr_i++)
    405                         {       size_t j = *itr_i;
    406                                 assert( j < f_.Domain() );
    407                                 std::set<size_t>::const_iterator itr_j;
     441                        entire_jac_sparse_.begin(i);
     442                        size_t j = entire_jac_sparse_.next_element();
     443                        while(j < n )
     444                        {       std::set<size_t>::const_iterator itr_j;
    408445                                const std::set<size_t>& r_j( r[j] );
    409446                                for(itr_j = r_j.begin(); itr_j != r_j.end(); itr_j++)
     
    412449                                        s[i].insert(k);
    413450                                }
     451                                j = entire_jac_sparse_.next_element();
    414452                        }
    415453                }
     
    430468                bool ok = true;
    431469                size_t m = f_.Range();
     470                size_t n = f_.Domain();
    432471                for(size_t i = 0; i < m; i++)
    433472                {       for(size_t k = 0; k < q; k++)
     
    438477                for(size_t i = 0; i < m; i++)
    439478                {       // compute row i of the return pattern
    440                         std::set<size_t>::const_iterator itr_i;
    441                         const std::set<size_t>& jac_i( entire_jac_sparse_[i] );
    442                         for(itr_i = jac_i.begin(); itr_i != jac_i.end(); itr_i++)
    443                         {       size_t j = *itr_i;
    444                                 assert( j < f_.Domain() );
    445                                 for(size_t k = 0; k < q; k++)
     479                        entire_jac_sparse_.begin(i);
     480                        size_t j = entire_jac_sparse_.next_element();
     481                        while(j < n )
     482                        {       for(size_t k = 0; k < q; k++)
    446483                                        s[i * q + k] |= r[j * q + k ];
     484                                j = entire_jac_sparse_.next_element();
    447485                        }
    448486                }
     
    463501                bool ok  = true;
    464502                //
     503                size_t m = f_.Range();
    465504                size_t n = f_.Domain();
    466505                for(size_t j = 0; j < n; j++)
     
    469508                // sparsity for  s = r * entire_jac_sparse_
    470509                // s^T = entire_jac_sparse_^T * r^T
    471                 for(size_t k = 0; k < q; k++)
    472                 {       // compute row k of the return pattern s
    473                         std::set<size_t>::const_iterator itr_k;
    474                         const std::set<size_t>& r_k( rt[k] );
    475                         for(itr_k = r_k.begin(); itr_k != r_k.end(); itr_k++)
    476                         {       size_t i = *itr_k;
    477                                 assert( i < f_.Range() );
    478                                 std::set<size_t>::const_iterator itr_i;
    479                                 const std::set<size_t>& jac_i( entire_jac_sparse_[i] );
    480                                 for(itr_i = jac_i.begin(); itr_i != jac_i.end(); itr_i++)
    481                                 {       size_t j = *itr_i;
    482                                         assert( j < n );
     510                for(size_t i = 0; i < m; i++)
     511                {       // i is the row index in r^T
     512                        std::set<size_t>::const_iterator itr_i;
     513                        const std::set<size_t>& r_i( rt[i] );
     514                        for(itr_i = r_i.begin(); itr_i != r_i.end(); itr_i++)
     515                        {       // k is the column index in r^T
     516                                size_t k = *itr_i;
     517                                assert( k < q );
     518                                //
     519                                // i is column index in entire_sparse_jac^T
     520                                entire_jac_sparse_.begin(i);
     521                                size_t j = entire_jac_sparse_.next_element();
     522                                while( j < n )
     523                                {       // j is row index in entire_sparse_jac^T
    483524                                        st[j].insert(k);
     525                                        j = entire_jac_sparse_.next_element();
    484526                                }
    485527                        }
     
    514556                        for(size_t i = 0; i < m; i++)
    515557                        {       if( rt[i * q + k] )
    516                                 {       std::set<size_t>::const_iterator itr_i;
    517                                         const std::set<size_t>& jac_i( entire_jac_sparse_[i] );
    518                                         for(itr_i = jac_i.begin(); itr_i != jac_i.end(); itr_i++)
    519                                         {       size_t j = *itr_i;
    520                                                 assert( j < n );
    521                                                 st[j * q + k ] = true;
     558                                {       entire_jac_sparse_.begin(i);
     559                                        size_t j = entire_jac_sparse_.next_element();
     560                                        while( j < n )
     561                                        {       st[j * q + k ] = true;
     562                                                j = entire_jac_sparse_.next_element();
    522563                                        }
    523564                                }
     
    551592                bool ok        = true;
    552593
    553                 // make sure entire_hes_sparse_ is set
    554                 if( entire_hes_sparse_.size() == 0 )
    555                 {
    556                         // set version of sparsity for vector of all ones
    557                         vector< std::set<size_t> > all_one(1);
    558                         CPPAD_ASSERT_UNKNOWN( all_one[0].empty() );
    559                         for(size_t i = 0; i < m; i++)
    560                                 all_one[0].insert(i);
    561 
    562                         // set version of sparsity for n by n idendity matrix
    563                         vector< std::set<size_t> > identity(n);
    564                         for(size_t j = 0; j < n; j++)
    565                                 identity[j].insert(j);
    566 
    567                         // compute sparsity pattern for H(x) = sum_i f_i(x)^{(2)}
    568                         bool transpose  = false;
    569                         bool dependency = false;
    570                         f_.ForSparseJac(n, identity, transpose, dependency);
    571                         entire_hes_sparse_ = f_.RevSparseHes(n, all_one, transpose);
    572 
    573                         // drop the forward sparsity results from f_
    574                         f_.size_forward_set(0);
    575                 }
     594                // make sure entire_hes_sparse_ has been set
     595                if( entire_hes_sparse_.n_set() == 0 )
     596                        set_entire_hes_sparse();
    576597
    577598                // compute sparsity pattern for T(x) = S(x) * f'(x)
     
    596617                for(size_t i = 0; i < n; i++)
    597618                {       v[i].clear();
    598                         std::set<size_t>::const_iterator itr_i;
    599                         const std::set<size_t>& hes_i( entire_hes_sparse_[i] );
    600                         for(itr_i = hes_i.begin(); itr_i != hes_i.end(); itr_i++)
    601                         {       size_t j = *itr_i;
    602                                 assert( j < n );
    603                                 std::set<size_t>::const_iterator itr_j;
     619                        entire_hes_sparse_.begin(i);
     620                        size_t j = entire_hes_sparse_.next_element();
     621                        while( j < n )
     622                        {       std::set<size_t>::const_iterator itr_j;
    604623                                const std::set<size_t>& r_j( r[j] );
    605624                                for(itr_j = r_j.begin(); itr_j != r_j.end(); itr_j++)
     
    607626                                        v[i].insert(k);
    608627                                }
     628                                j = entire_hes_sparse_.next_element();
    609629                        }
    610630                }
     
    645665                bool ok        = true;
    646666
    647                 // make sure entire_hes_sparse_ is set
    648                 if( entire_hes_sparse_.size() == 0 )
    649                 {
    650                         // set version of sparsity for vector of all ones
    651                         vector< std::set<size_t> > all_one(1);
    652                         CPPAD_ASSERT_UNKNOWN( all_one[0].empty() );
    653                         for(size_t i = 0; i < m; i++)
    654                                 all_one[0].insert(i);
    655 
    656                         // set version of sparsity for n by n idendity matrix
    657                         vector< std::set<size_t> > identity(n);
    658                         for(size_t j = 0; j < n; j++)
    659                                 identity[j].insert(j);
    660 
    661                         // compute sparsity pattern for H(x) = sum_i f_i(x)^{(2)}
    662                         bool transpose  = false;
    663                         bool dependency = false;
    664                         f_.ForSparseJac(n, identity, transpose, dependency);
    665                         entire_hes_sparse_ = f_.RevSparseHes(n, all_one, transpose);
    666 
    667                         // drop the forward sparsity results from f_
    668                         f_.size_forward_set(0);
    669                 }
     667                // make sure entire_hes_sparse_ has been set
     668                if( entire_hes_sparse_.n_set() == 0 )
     669                        set_entire_hes_sparse();
    670670
    671671                // compute sparsity pattern for T(x) = S(x) * f'(x)
     
    691691                {       for(size_t k = 0; k < q; k++)
    692692                                v[i * q + k] = false;
    693                         std::set<size_t>::const_iterator itr_i;
    694                         const std::set<size_t>& hes_i( entire_hes_sparse_[i] );
    695                         for(itr_i = hes_i.begin(); itr_i != hes_i.end(); itr_i++)
    696                         {       size_t j = *itr_i;
    697                                 assert( j < n );
    698                                 for(size_t k = 0; k < q; k++)
     693                        entire_hes_sparse_.begin(i);
     694                        size_t j = entire_hes_sparse_.next_element();
     695                        while( j < n )
     696                        {       for(size_t k = 0; k < q; k++)
    699697                                        v[i * q + k] |= r[ j * q + k ];
     698                                j = entire_hes_sparse_.next_element();
    700699                        }
    701700                }
  • trunk/cppad/local/for_sparse_jac.hpp

    r3711 r3712  
    755755        return s;
    756756}
     757// ===========================================================================
     758// ForSparseJacCheckpoint
     759/*!
     760Forward mode Jacobian sparsity calculation used by checkpoint functions.
     761
     762\tparam Base
     763is the base type for this recording.
     764
     765\param transpose
     766is true (false) s is equal to \f$ S(x) \f$ (\f$ S(x)^T \f$)
     767where
     768\f[
     769        S(x) = F^{(1)} (x) * R
     770\f]
     771where \f$ F \f$ is the function corresponding to the operation sequence
     772and \f$ x \f$ is any argument value.
     773
     774\param q
     775is the number of columns in the matrix \f$ R \f$.
     776
     777\param r
     778is a sparsity pattern for the matrix \f$ R \f$.
     779
     780\param transpose
     781are the sparsity patterns for \f$ R \f$ and \f$ S(x) \f$ transposed.
     782
     783\param dependency
     784Are the derivatives with respect to left and right of the expression below
     785considered to be non-zero:
     786\code
     787        CondExpRel(left, right, if_true, if_false)
     788\endcode
     789This is used by the optimizer to obtain the correct dependency relations.
     790
     791\param s
     792The input size and elements of s do not matter.
     793On output, s is the sparsity pattern for the matrix \f$ S(x) \f$
     794or \f$ S(x)^T \f$ depending on transpose.
     795*/
     796template <class Base>
     797void ADFun<Base>::ForSparseJacCheckpoint(
     798        size_t                        q          ,
     799        CPPAD_INTERNAL_SPARSE_SET&    r          ,
     800        bool                          transpose  ,
     801        bool                          dependency ,
     802        CPPAD_INTERNAL_SPARSE_SET&    s          )
     803{       size_t n = Domain();
     804        size_t m = Range();
     805
     806# ifndef NDEBUG
     807        if( transpose )
     808        {       CPPAD_ASSERT_UNKNOWN( r.n_set() == q );
     809                CPPAD_ASSERT_UNKNOWN( r.end()   == n );
     810        }
     811        else
     812        {       CPPAD_ASSERT_UNKNOWN( r.n_set() == n );
     813                CPPAD_ASSERT_UNKNOWN( r.end()   == q );
     814        }
     815        for(size_t j = 0; j < n; j++)
     816        {       CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] == (j+1) );
     817                CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[j] ) == InvOp );
     818        }
     819# endif
     820
     821        // holds reverse Jacobian sparsity pattern for all variables
     822        CPPAD_INTERNAL_SPARSE_SET var_sparsity;
     823        var_sparsity.resize(num_var_tape_, q);
     824
     825        // set sparsity pattern for dependent variables
     826        if( transpose )
     827        {       for(size_t i = 0; i < q; i++)
     828                {       r.begin(i);
     829                        size_t j = r.next_element();
     830                        while( j < n )
     831                        {       var_sparsity.add_element( ind_taddr_[j], i );
     832                                j = r.next_element();
     833                        }
     834                }
     835        }
     836        else
     837        {       for(size_t j = 0; j < n; j++)
     838                {       r.begin(j);
     839                        size_t i = r.next_element();
     840                        while( i < q )
     841                        {       var_sparsity.add_element( ind_taddr_[j], i );
     842                                i = r.next_element();
     843                        }
     844                }
     845        }
     846
     847        // evaluate the sparsity pattern for all variables
     848        ForJacSweep(
     849                dependency,
     850                n,
     851                num_var_tape_,
     852                &play_,
     853                var_sparsity
     854        );
     855
     856        // dimension the return value
     857        if( transpose )
     858                s.resize(q, m);
     859        else
     860                s.resize(m, q);
     861
     862        // return values corresponding to dependent variables
     863        for(size_t i = 0; i < m; i++)
     864        {       CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < num_var_tape_ );
     865
     866                // extract the result from var_sparsity
     867                CPPAD_ASSERT_UNKNOWN( var_sparsity.end() == q );
     868                var_sparsity.begin( dep_taddr_[i] );
     869                size_t j = var_sparsity.next_element();
     870                while( j < q )
     871                {       if( transpose )
     872                                s.add_element(j, i);
     873                        else
     874                                s.add_element(i, j);
     875                        j  = var_sparsity.next_element();
     876                }
     877        }
     878
     879}
    757880
    758881
  • trunk/cppad/local/rev_sparse_hes.hpp

    r3711 r3712  
    200200*/
    201201
     202// ===========================================================================
     203// RevSparseHesBool
    202204/*!
    203205Calculate Hessian sparsity patterns using reverse mode.
     
    355357}
    356358
     359// ===========================================================================
     360// RevSparseHesSet
    357361/*!
    358362Calculate Hessian sparsity patterns using reverse mode.
     
    515519}
    516520
    517 
     521// ===========================================================================
     522// RevSparseHes
    518523
    519524/*!
     
    587592        return h;
    588593}
    589 
     594// ===========================================================================
     595// RevSparseHesCase
    590596/*!
    591597Private helper function for RevSparseHes(q, s).
     
    642648        );
    643649}
    644 
    645650/*!
    646651Private helper function for RevSparseHes(q, s).
     
    700705        );
    701706}
     707// ===========================================================================
     708// RevSparseHesCheckpoint
     709/*!
     710Hessian sparsity patterns calculation used by checkpoint functions.
     711
     712\tparam Base
     713is the base type for this recording.
     714
     715\param transpose
     716is true (false) h is equal to \f$ H(x) \f$ (\f$ H(x)^T \f$)
     717where
     718\f[
     719        H(x) = R^T (S * F)^{(2)} (x)
     720\f]
     721where \f$ F \f$ is the function corresponding to the operation sequence
     722and \f$ x \f$ is any argument value.
     723
     724\param q
     725is the value of q in the by the previous call of the form
     726\verbatim
     727        f.ForSparseJac(q, r)
     728\endverbatim
     729The value r in this call is a sparsity pattern for the matrix \f$ R \f$.
     730
     731\param s
     732is a vector with size m that specifies the sparsity pattern
     733for the vector \f$ S \f$,
     734where m is the number of dependent variables
     735corresponding to the operation sequence stored in play_.
     736
     737\param h
     738The input size and elements of h do not matter.
     739On output, h is the sparsity pattern for the matrix \f$ H(x) \f$
     740or \f$ H(x)^T \f$ depending on transpose.
     741
     742\par Assumptions
     743The forward jacobian sparsity pattern must be currently stored
     744in this ADFUN object.
     745*/
     746template <class Base>
     747void ADFun<Base>::RevSparseHesCheckpoint(
     748        size_t                        q         ,
     749        vector<bool>&                 s         ,
     750        bool                          transpose ,
     751        CPPAD_INTERNAL_SPARSE_SET&    h         )
     752{       size_t n = Domain();
     753        size_t m = Range();
     754
     755        // checkpoint functions should get this right
     756        CPPAD_ASSERT_UNKNOWN( for_jac_sparse_pack_.n_set() == 0 );
     757        CPPAD_ASSERT_UNKNOWN( for_jac_sparse_set_.n_set() == num_var_tape_ );
     758        CPPAD_ASSERT_UNKNOWN( for_jac_sparse_set_.end()   == q );
     759        CPPAD_ASSERT_UNKNOWN( s.size()                    == m );
     760
     761        // Array that holds the reverse Jacobiain dependcy flags.
     762        // Initialize as true for dependent variables, flase for others.
     763        pod_vector<bool> RevJac;
     764        RevJac.extend(num_var_tape_);
     765        for(size_t i = 0; i < num_var_tape_; i++)
     766                RevJac[i] = false;
     767        for(size_t i = 0; i < m; i++)
     768        {       CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < num_var_tape_ )
     769                RevJac[ dep_taddr_[i] ] = s[i];
     770        }
     771
     772        // holds reverse Hessian sparsity pattern for all variables
     773        CPPAD_INTERNAL_SPARSE_SET rev_hes_sparsity;
     774        rev_hes_sparsity.resize(num_var_tape_, q);
     775
     776        // compute Hessian sparsity pattern for all variables
     777        RevHesSweep(
     778                n,
     779                num_var_tape_,
     780                &play_,
     781                for_jac_sparse_set_,
     782                RevJac.data(),
     783                rev_hes_sparsity
     784        );
     785
     786        // dimension the return value
     787        if( transpose )
     788                h.resize(n, q);
     789        else
     790                h.resize(q, n);
     791
     792        // j is index corresponding to reverse mode partial
     793        for(size_t j = 0; j < n; j++)
     794        {       CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] < num_var_tape_ );
     795
     796                // ind_taddr_[j] is operator taddr for j-th independent variable
     797                CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] == j + 1 );
     798                CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[j] ) == InvOp );
     799
     800                // extract the result from rev_hes_sparsity
     801                CPPAD_ASSERT_UNKNOWN( rev_hes_sparsity.end() == q );
     802                rev_hes_sparsity.begin(j + 1);
     803                size_t i = rev_hes_sparsity.next_element();
     804                while( i < q )
     805                {       if( transpose )
     806                                h.add_element(j,  i);
     807                        else    h.add_element(i, j);
     808                        i = rev_hes_sparsity.next_element();
     809                }
     810        }
     811}
    702812
    703813} // END_CPPAD_NAMESPACE
  • trunk/cppad/local/rev_sparse_jac.hpp

    r3710 r3712  
    191191*/
    192192
    193 // -------------------------------------------------------------------------
     193// =========================================================================
     194// RevSparseJacBool
    194195/*!
    195196Calculate Jacobian vector of bools sparsity patterns using reverse mode.
     
    337338        }
    338339}
     340// =========================================================================
     341// RevSparseJacSet
    339342/*!
    340343Calculate Jacobian vector of sets sparsity patterns using reverse mode.
     
    480483        }
    481484}
    482 // --------------------------------------------------------------------------
     485// =========================================================================
     486// RevSparseJacCase
    483487
    484488/*!
     
    592596        );
    593597}
    594 // --------------------------------------------------------------------------
     598// =========================================================================
     599// RevSparseJac
    595600/*!
    596601User API for Jacobian sparsity patterns using reverse mode.
     
    661666        return s;
    662667}
     668// ===========================================================================
     669// RevSparseJacCheckpoint
     670/*!
     671Reverse mode Jacobian sparsity calculation used by checkpoint functions.
     672
     673\tparam Base
     674is the base type for this recording.
     675
     676\param transpose
     677is true (false) s is equal to \f$ S(x) \f$ (\f$ S(x)^T \f$)
     678where
     679\f[
     680        S(x) = R * F^{(1)} (x)
     681\f]
     682where \f$ F \f$ is the function corresponding to the operation sequence
     683and \f$ x \f$ is any argument value.
     684
     685\param q
     686is the number of rows in the matrix \f$ R \f$.
     687
     688\param r
     689is a sparsity pattern for the matrix \f$ R \f$.
     690
     691\param transpose
     692are the sparsity patterns for \f$ R \f$ and \f$ S(x) \f$ transposed.
     693
     694\param dependency
     695Are the derivatives with respect to left and right of the expression below
     696considered to be non-zero:
     697\code
     698        CondExpRel(left, right, if_true, if_false)
     699\endcode
     700This is used by the optimizer to obtain the correct dependency relations.
     701
     702\param s
     703The input size and elements of s do not matter.
     704On output, s is the sparsity pattern for the matrix \f$ S(x) \f$
     705or \f$ S(x)^T \f$ depending on transpose.
     706
     707*/
     708template <class Base>
     709void ADFun<Base>::RevSparseJacCheckpoint(
     710        size_t                        q          ,
     711        CPPAD_INTERNAL_SPARSE_SET&    r          ,
     712        bool                          transpose  ,
     713        bool                          dependency ,
     714        CPPAD_INTERNAL_SPARSE_SET&    s          )
     715{       size_t n = Domain();
     716        size_t m = Range();
     717
     718# ifndef NDEBUG
     719        if( transpose )
     720        {       CPPAD_ASSERT_UNKNOWN( r.n_set() == m );
     721                CPPAD_ASSERT_UNKNOWN( r.end()   == q );
     722        }
     723        else
     724        {       CPPAD_ASSERT_UNKNOWN( r.n_set() == q );
     725                CPPAD_ASSERT_UNKNOWN( r.end()   == m );
     726        }
     727        for(size_t i = 0; i < m; i++)
     728                CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < num_var_tape_ );
     729# endif
     730
     731        // holds reverse Jacobian sparsity pattern for all variables
     732        CPPAD_INTERNAL_SPARSE_SET var_sparsity;
     733        var_sparsity.resize(num_var_tape_, q);
     734
     735        // set sparsity pattern for dependent variables
     736        if( transpose )
     737        {       for(size_t i = 0; i < m; i++)
     738                {       r.begin(i);
     739                        size_t j = r.next_element();
     740                        while( j < q )
     741                        {       var_sparsity.add_element( dep_taddr_[i], j );
     742                                j = r.next_element();
     743                        }
     744                }
     745        }
     746        else
     747        {       for(size_t j = 0; j < q; j++)
     748                {       r.begin(j);
     749                        size_t i = r.next_element();
     750                        while( i < m )
     751                        {       var_sparsity.add_element( dep_taddr_[i], j );
     752                                i = r.next_element();
     753                        }
     754                }
     755        }
     756
     757        // evaluate the sparsity pattern for all variables
     758        RevJacSweep(
     759                dependency,
     760                n,
     761                num_var_tape_,
     762                &play_,
     763                var_sparsity
     764        );
     765
     766        // dimension the return value
     767        if( transpose )
     768                s.resize(n, m);
     769        else
     770                s.resize(m, n);
     771
     772        // return values corresponding to independent variables
     773        for(size_t j = 0; j < n; j++)
     774        {       CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] == (j+1) );
     775
     776                // ind_taddr_[j] is operator taddr for j-th independent variable
     777                CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[j] ) == InvOp );
     778
     779                // extract the result from var_sparsity
     780                CPPAD_ASSERT_UNKNOWN( var_sparsity.end() == q );
     781                var_sparsity.begin(j+1);
     782                size_t i = var_sparsity.next_element();
     783                while( i < q )
     784                {       if( transpose )
     785                                s.add_element(j, i);
     786                        else
     787                                s.add_element(i, j);
     788                        i  = var_sparsity.next_element();
     789                }
     790        }
     791
     792}
    663793
    664794} // END_CPPAD_NAMESPACE
  • trunk/omh/whats_new/whats_new_15.omh

    r3711 r3712  
    6161assist you in learning about changes between various versions of CppAD.
    6262
     63$head 08-26$$
     64Fix a bug in $cref RevSparseJac$$ when used to compute sparsity pattern
     65for a subset of the rows in a $cref checkpoint$$ function.
     66
     67$head 08-25$$
     68Reduce the amount of memory required for $cref checkpoint$$ functions
     69(since sparsity patterns are now being held so they do not need to be
     70recalculated).
     71
    6372$head 08-20$$
    6473Added an example that computes the sparsity pattern for a subset
  • trunk/test_more/checkpoint.cpp

    r3708 r3712  
    11/* $Id$ */
    22/* --------------------------------------------------------------------------
    3 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
     3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-15 Bradley M. Bell
    44
    55CppAD is distributed under multiple licenses. This distribution is under
     
    135135                return ok;
    136136        }
     137
     138        bool h_algo(const ADVector& ax, ADVector& ay)
     139        {       ay[0] = ax[0];
     140                ay[1] = ax[1] + ax[2];
     141                return true;
     142        }
     143        bool test_two(void)
     144        {       bool ok = true;
     145                using CppAD::checkpoint;
     146                using CppAD::ADFun;
     147                using CppAD::NearEqual;
     148
     149                // checkpoint version of H(x)
     150                size_t m = 2;
     151                size_t n = 3;
     152                ADVector ax(n), ay(m);
     153                for(size_t j = 0; j < n; j++)
     154                        ax[j] = double(j);
     155                checkpoint<double> h_check("h_check", h_algo, ax, ay);
     156
     157                // record function using h_check
     158                Independent(ax);
     159                h_check(ax, ay);
     160                ADFun<double> h(ax, ay);
     161
     162                // compute sparsity pattern h_1(x) = x[1] + x[2]
     163                CppAD::vector< std::set<size_t> > r(1), s(1);
     164                r[0].insert(1);
     165                s = h.RevSparseJac(1, r);
     166
     167                // check result
     168                ok &= s[0] == std::set<size_t>{1, 2};
     169
     170                return ok;
     171        }
    137172}
    138173
     
    140175{       bool ok = true;
    141176        ok  &= test_one();
    142 
     177        ok  &= test_two();
    143178        return ok;
    144179}
Note: See TracChangeset for help on using the changeset viewer.