Ignore:
Timestamp:
Mar 25, 2016 12:36:53 AM (4 years ago)
Author:
bradbell
Message:

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

commit dd967ef41b8d6731d90ebb8b3e7e8b863565289c
Author: Brad Bell <bradbell@…>
Date: Thu Mar 24 19:19:00 2016 -0700

eigen_mat_mul.hpp: add for_sparse_hes calculation.
atomic_base.hpp: edits to forward and reverse sparse hessian documentation.
for_hes_sweep.hpp: fix bug in ForSparseHes? calculation.
eigen_mat_mul.cpp: test for_sparse_hes calculation.

commit 99610d3bdda1162ec7e8884b3cc5811b5fc48c24
Author: Brad Bell <bradbell@…>
Date: Thu Mar 24 17:51:42 2016 -0700

eigen_mat_mul.cpp: test rev_sparse_hes.
eigen_mat_mul.hpp: fix heading -> subheading.

commit 19a0f5d8c3c1e8ea210f852f49adc290609fdf10
Author: Brad Bell <bradbell@…>
Date: Thu Mar 24 17:24:50 2016 -0700

eigen_mat_mul.hpp: add code for rev_sparse_hes.
atomic_base.hpp: in doc change some g(y) -> g[f(x)] (clearer).
eigen_mat_mul.cpp: test second order derivatives.

commit 4487cc1d4f5e598d690dba681b5c275281981bf0
Author: Brad Bell <bradbell@…>
Date: Thu Mar 24 15:14:17 2016 -0700

Add rev_sparse_jac to eigen_mat_mul.hpp.

commit 055fa95218ca47e30204796c887f33b5ca9f9788
Author: Brad Bell <bradbell@…>
Date: Thu Mar 24 07:40:49 2016 -0700

eigen_mat_mul.hpp: use subheadings to separate Public and Private.
eigen_mat_mul.cpp: test for_sparse_jac.

commit af8898d93493df4d1a088038abe926f9ab9f54d5
Author: Brad Bell <bradbell@…>
Date: Thu Mar 24 07:06:46 2016 -0700

eigen_mat_mul.hpp: add for_sparse_jac (not yet tested).
eigen_mat_mul.cpp: change to example with a non-zero Hessian.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/cppad/example/eigen_mat_mul.hpp

    r3808 r3809  
    3131
    3232/* %$$
    33 $head Publice Types$$
     33$head Publice$$
     34
     35$subhead Types$$
    3436$srccode%cpp% */
    3537namespace { // BEGIN_EMPTY_NAMESPACE
     
    5052                ad_scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > ad_matrix;
    5153/* %$$
    52 $head Public Constructor$$
     54$subhead Constructor$$
    5355$srccode%cpp% */
    5456        // constructor
     
    7274        { }
    7375/* %$$
    74 $head Public Pack$$
     76$subhead Pack$$
    7577$srccode%cpp% */
    7678        template <class Matrix, class Vector>
     
    9799        }
    98100/* %$$
    99 $head Public Unpack$$
     101$subhead Unpack$$
    100102$srccode%cpp% */
    101103        template <class Matrix, class Vector>
     
    116118        }
    117119/* %$$
    118 $head Private Variables$$
     120$head Private$$
     121
     122$subhead Variables$$
    119123$srccode%cpp% */
    120124private:
     
    137141        // -------------------------------------------------------------
    138142/* %$$
    139 $head Private rows$$
     143$subhead rows$$
    140144$srccode%cpp% */
    141145        // convert from int to size_t
     
    145149        {       return size_t( x.rows() ); }
    146150/* %$$
    147 $head Private cols$$
     151$subhead cols$$
    148152$srccode%cpp% */
    149153        // convert from int to size_t
     
    153157        {       return size_t( x.cols() ); }
    154158/* %$$
    155 $head Private forward$$
     159$subhead forward$$
    156160$srccode%cpp% */
    157161        // forward mode routine called by CppAD
     
    252256        }
    253257/* %$$
    254 $head Private reverse$$
     258$subhead reverse$$
    255259$srccode%cpp% */
    256260        // reverse mode routine called by CppAD
     
    343347        }
    344348/* %$$
     349$subhead for_sparse_jac$$
     350$srccode%cpp% */
     351        // forward Jacobian sparsity routine called by CppAD
     352        virtual bool for_sparse_jac(
     353                // number of columns in the matrix R
     354                size_t                                       q ,
     355                // sparsity pattern for the matrix R
     356                const CppAD::vector< std::set<size_t> >&     r ,
     357                // sparsity pattern for the matrix S = f'(x) * R
     358                CppAD::vector< std::set<size_t> >&           s )
     359        {       assert( nx_ == r.size() );
     360                assert( ny_ == s.size() );
     361                //
     362                size_t n_left = nr_left_ * n_middle_;
     363                for(size_t i = 0; i < nr_left_; i++)
     364                {       for(size_t j = 0; j < nc_right_; j++)
     365                        {       // pack index for entry (i, j) in result
     366                                size_t i_result = i * nc_right_ + j;
     367                                s[i_result].clear();
     368                                for(size_t ell = 0; ell < n_middle_; ell++)
     369                                {       // pack index for entry (i, ell) in left
     370                                        size_t i_left  = i * n_middle_ + ell;
     371                                        // pack index for entry (ell, j) in right
     372                                        size_t i_right = n_left + ell * nc_right_ + j;
     373                                        //
     374                                        s[i_result] = CppAD::set_union(s[i_result], r[i_left] );
     375                                        s[i_result] = CppAD::set_union(s[i_result], r[i_right] );
     376                                }
     377                        }
     378                }
     379                return true;
     380        }
     381/* %$$
     382$subhead rev_sparse_jac$$
     383$srccode%cpp% */
     384        // reverse Jacobian sparsity routine called by CppAD
     385        virtual bool rev_sparse_jac(
     386                // number of columns in the matrix R^T
     387                size_t                                      q  ,
     388                // sparsity pattern for the matrix R^T
     389                const CppAD::vector< std::set<size_t> >&    rt ,
     390                // sparsoity pattern for the matrix S^T = f'(x)^T * R^T
     391                CppAD::vector< std::set<size_t> >&          st )
     392        {       assert( nx_ == st.size() );
     393                assert( ny_ == rt.size() );
     394
     395                // initialize S^T as empty
     396                for(size_t i = 0; i < nx_; i++)
     397                        st[i].clear();
     398
     399                // sparsity for S(x)^T = f'(x)^T * R^T
     400                size_t n_left = nr_left_ * n_middle_;
     401                for(size_t i = 0; i < nr_left_; i++)
     402                {       for(size_t j = 0; j < nc_right_; j++)
     403                        {       // pack index for entry (i, j) in result
     404                                size_t i_result = i * nc_right_ + j;
     405                                st[i_result].clear();
     406                                for(size_t ell = 0; ell < n_middle_; ell++)
     407                                {       // pack index for entry (i, ell) in left
     408                                        size_t i_left  = i * n_middle_ + ell;
     409                                        // pack index for entry (ell, j) in right
     410                                        size_t i_right = n_left + ell * nc_right_ + j;
     411                                        //
     412                                        st[i_left]  = CppAD::set_union(st[i_left],  rt[i_result]);
     413                                        st[i_right] = CppAD::set_union(st[i_right], rt[i_result]);
     414                                }
     415                        }
     416                }
     417                return true;
     418        }
     419/* %$$
     420$subhead for_sparse_hes$$
     421$srccode%cpp% */
     422        virtual bool for_sparse_hes(
     423                // which components of x are variables for this call
     424                const CppAD::vector<bool>&                   vx,
     425                // sparsity pattern for the diagonal of R
     426                const CppAD::vector<bool>&                   r ,
     427                // sparsity pattern for the vector S
     428                const CppAD::vector<bool>&                   s ,
     429                // sparsity patternfor the Hessian H(x)
     430                CppAD::vector< std::set<size_t> >&           h )
     431        {       assert( vx.size() == nx_ );
     432                assert( r.size()  == nx_ );
     433                assert( s.size()  == ny_ );
     434                assert( h.size()  == nx_ );
     435                //
     436                // initilize h as empty
     437                for(size_t i = 0; i < nx_; i++)
     438                        h[i].clear();
     439                //
     440                size_t n_left = nr_left_ * n_middle_;
     441                for(size_t i = 0; i < nr_left_; i++)
     442                {       for(size_t j = 0; j < nc_right_; j++)
     443                        {       // pack index for entry (i, j) in result
     444                                size_t i_result = i * nc_right_ + j;
     445                                if( s[i_result] )
     446                                {       for(size_t ell = 0; ell < n_middle_; ell++)
     447                                        {       // pack index for entry (i, ell) in left
     448                                                size_t i_left  = i * n_middle_ + ell;
     449                                                // pack index for entry (ell, j) in right
     450                                                size_t i_right = n_left + ell * nc_right_ + j;
     451                                                if( r[i_left] & r[i_right] )
     452                                                {       h[i_left].insert(i_right);
     453                                                        h[i_right].insert(i_left);
     454                                                }
     455                                        }
     456                                }
     457                        }
     458                }
     459                return true;
     460        }
     461/* %$$
     462$subhead rev_sparse_hes$$
     463$srccode%cpp% */
     464        // reverse Hessian sparsity routine called by CppAD
     465        virtual bool rev_sparse_hes(
     466                // which components of x are variables for this call
     467                const CppAD::vector<bool>&                   vx,
     468                // sparsity pattern for S(x) = g'[f(x)]
     469                const CppAD::vector<bool>&                   s ,
     470                // sparsity pattern for d/dx g[f(x)] = S(x) * f'(x)
     471                CppAD::vector<bool>&                         t ,
     472                // number of columns in R, U(x), and V(x)
     473                size_t                                       q ,
     474                // sparsity pattern for R
     475                const CppAD::vector< std::set<size_t> >&     r ,
     476                // sparsity pattern for U(x) = g^{(2)} [ f(x) ] * f'(x) * R
     477                const CppAD::vector< std::set<size_t> >&     u ,
     478                // sparsity pattern for
     479                // V(x) = f'(x)^T * U(x) + sum_{i=0}^{m-1} S_i(x) f_i^{(2)} (x) * R
     480                CppAD::vector< std::set<size_t> >&           v )
     481        {       assert( vx.size() == nx_ );
     482                assert( s.size()  == ny_ );
     483                assert( t.size()  == nx_ );
     484                assert( r.size()  == nx_ );
     485                assert( v.size()  == nx_ );
     486                //
     487                // initilaize return sparsity patterns as false
     488                for(size_t j = 0; j < nx_; j++)
     489                {       t[j] = false;
     490                        v[j].clear();
     491                }
     492                //
     493                size_t n_left = nr_left_ * n_middle_;
     494                for(size_t i = 0; i < nr_left_; i++)
     495                {       for(size_t j = 0; j < nc_right_; j++)
     496                        {       // pack index for entry (i, j) in result
     497                                size_t i_result = i * nc_right_ + j;
     498                                for(size_t ell = 0; ell < n_middle_; ell++)
     499                                {       // pack index for entry (i, ell) in left
     500                                        size_t i_left  = i * n_middle_ + ell;
     501                                        // pack index for entry (ell, j) in right
     502                                        size_t i_right = n_left + ell * nc_right_ + j;
     503                                        //
     504                                        // back propagate T(x) = S(x) * f'(x).
     505                                        t[i_left]  |= bool( s[i_result] );
     506                                        t[i_right] |= bool( s[i_result] );
     507                                        //
     508                                        // V(x) = f'(x)^T * U(x) +  sum_i S_i(x) * f_i''(x) * R
     509                                        // U(x)   = g''[ f(x) ] * f'(x) * R
     510                                        // S_i(x) = g_i'[ f(x) ]
     511                                        //
     512                                        // back propagate f'(x)^T * U(x)
     513                                        v[i_left]  = CppAD::set_union(v[i_left],  u[i_result] );
     514                                        v[i_right] = CppAD::set_union(v[i_right], u[i_result] );
     515                                        //
     516                                        // back propagate S_i(x) * f_i''(x) * R
     517                                        // (here is where we use vx to check for cross terms)
     518                                        if( s[i_result] & vx[i_left] & vx[i_right] )
     519                                        {       v[i_left]  = CppAD::set_union(v[i_left],  r[i_right] );
     520                                                v[i_right] = CppAD::set_union(v[i_right], r[i_left]  );
     521                                        }
     522                                }
     523                        }
     524                }
     525                return true;
     526        }
     527/* %$$
    345528$head End Class Definition$$
    346529$srccode%cpp% */
Note: See TracChangeset for help on using the changeset viewer.