source: trunk/cppad/example/eigen_mat_mul.hpp @ 3809

Last change on this file since 3809 was 3809, checked in by bradbell, 4 years ago

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 size: 17.2 KB
Line 
1// $Id$
2# ifndef CPPAD_EXAMPLE_EIGEN_MAT_MUL_HPP
3# define CPPAD_EXAMPLE_EIGEN_MAT_MUL_HPP
4
5/* --------------------------------------------------------------------------
6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-16 Bradley M. Bell
7
8CppAD is distributed under multiple licenses. This distribution is under
9the terms of the
10                    Eclipse Public License Version 1.0.
11
12A copy of this license is included in the COPYING file of this distribution.
13Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
14-------------------------------------------------------------------------- */
15
16/*
17$begin atomic_eigen_mat_mul.hpp$$
18$spell
19        Eigen
20$$
21
22$section Atomic Eigen Matrix Multiply Class$$
23
24$nospell
25
26$head Start Class Definition$$
27$srccode%cpp% */
28# include <cppad/cppad.hpp>
29# include <Eigen/Core>
30
31
32/* %$$
33$head Publice$$
34
35$subhead Types$$
36$srccode%cpp% */
37namespace { // BEGIN_EMPTY_NAMESPACE
38
39template <class Base>
40class atomic_eigen_mat_mul : public CppAD::atomic_base<Base> {
41public:
42        // -----------------------------------------------------------
43        // type of elements during calculation of derivatives
44        typedef Base              scalar;
45        // type of elements during taping
46        typedef CppAD::AD<scalar> ad_scalar;
47        // type of matrix during calculation of derivatives
48        typedef Eigen::Matrix<
49                scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>     matrix;
50        // type of matrix during taping
51        typedef Eigen::Matrix<
52                ad_scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > ad_matrix;
53/* %$$
54$subhead Constructor$$
55$srccode%cpp% */
56        // constructor
57        atomic_eigen_mat_mul(
58                // number of rows in left operand
59                const size_t nr_left  ,
60                // number of rows in left and columns in right operand
61                const size_t n_middle   ,
62                // number of columns in right operand
63                const size_t nc_right
64        ) :
65        CppAD::atomic_base<Base>(
66                "atom_eigen_mat_mul"                             ,
67                CppAD::atomic_base<Base>::set_sparsity_enum
68        ) ,
69        nr_left_(  nr_left  )                   ,
70        n_middle_(   n_middle   )                   ,
71        nc_right_( nc_right )                   ,
72        nx_( (nr_left + nc_right) * n_middle_ ) ,
73        ny_( nr_left * nc_right )
74        { }
75/* %$$
76$subhead Pack$$
77$srccode%cpp% */
78        template <class Matrix, class Vector>
79        void pack(
80                Vector&        packed  ,
81                const Matrix&  left    ,
82                const Matrix&  right   )
83        {       assert( packed.size() == nx_      );
84                assert( rows( left ) == nr_left_ );
85                assert( cols( left ) == n_middle_ );
86                assert( rows( right ) == n_middle_ );
87                assert( cols( right ) == nc_right_ );
88                //
89                size_t n_left = nr_left_ * n_middle_;
90                size_t n_right = n_middle_ * nc_right_;
91                assert( n_left + n_right == nx_ );
92                //
93                for(size_t i = 0; i < n_left; i++)
94                        packed[i] = left.data()[i];
95                for(size_t i = 0; i < n_right; i++)
96                        packed[ i + n_left ] = right.data()[i];
97                //
98                return;
99        }
100/* %$$
101$subhead Unpack$$
102$srccode%cpp% */
103        template <class Matrix, class Vector>
104        void unpack(
105                const Vector&   packed  ,
106                Matrix&         result  )
107        {       assert( packed.size() == ny_      );
108                assert( rows( result ) == nr_left_ );
109                assert( cols( result ) == nc_right_ );
110                //
111                size_t n_result = nr_left_ * nc_right_;
112                assert( n_result == ny_ );
113                //
114                for(size_t i = 0; i < n_result; i++)
115                        result.data()[i] = packed[ i ];
116                //
117                return;
118        }
119/* %$$
120$head Private$$
121
122$subhead Variables$$
123$srccode%cpp% */
124private:
125        // -------------------------------------------------------------
126        // number of of rows in left operand and result
127        const size_t nr_left_;
128        // number of of columns in left operand and rows in right operand
129        const size_t n_middle_;
130        // number of columns in right operand and result
131        const size_t nc_right_;
132        // size of the domain space
133        const size_t nx_;
134        // size of the range space
135        const size_t ny_;
136        // -------------------------------------------------------------
137        // one forward mode vector of matrices for left, right, and result
138        CppAD::vector<matrix> f_left_, f_right_, f_result_;
139        // one reverse mode vector of matrices for left, right, and result
140        CppAD::vector<matrix> r_left_, r_right_, r_result_;
141        // -------------------------------------------------------------
142/* %$$
143$subhead rows$$
144$srccode%cpp% */
145        // convert from int to size_t
146        static size_t rows(const matrix& x)
147        {       return size_t( x.rows() ); }
148        static size_t rows(const ad_matrix& x)
149        {       return size_t( x.rows() ); }
150/* %$$
151$subhead cols$$
152$srccode%cpp% */
153        // convert from int to size_t
154        static size_t cols(const matrix& x)
155        {       return size_t( x.cols() ); }
156        static size_t cols(const ad_matrix& x)
157        {       return size_t( x.cols() ); }
158/* %$$
159$subhead forward$$
160$srccode%cpp% */
161        // forward mode routine called by CppAD
162        virtual bool forward(
163                // lowest order Taylor coefficient we are evaluating
164                size_t                          p ,
165                // highest order Taylor coefficient we are evaluating
166                size_t                          q ,
167                // which components of x are variables
168                const CppAD::vector<bool>&      vx ,
169                // which components of y are variables
170                CppAD::vector<bool>&            vy ,
171                // tx [ j * (q+1) + k ] is x_j^k
172                const CppAD::vector<scalar>&    tx ,
173                // ty [ i * (q+1) + k ] is y_i^k
174                CppAD::vector<scalar>&          ty
175        )
176        {       size_t n_order = q + 1;
177                assert( vx.size() == 0 || nx_ == vx.size() );
178                assert( vx.size() == 0 || ny_ == vy.size() );
179                assert( nx_ * n_order == tx.size() );
180                assert( ny_ * n_order == ty.size() );
181                //
182                size_t n_left   = nr_left_ * n_middle_;
183                size_t n_right  = n_middle_ * nc_right_;
184                size_t n_result = nr_left_ * nc_right_;
185                assert( n_left + n_right == nx_ );
186                assert( n_result == ny_ );
187                //
188                // -------------------------------------------------------------------
189                // make sure f_left_, f_right_, and f_result_ are large enough
190                assert( f_left_.size() == f_right_.size() );
191                assert( f_left_.size() == f_result_.size() );
192                if( f_left_.size() < n_order )
193                {       f_left_.resize(n_order);
194                        f_right_.resize(n_order);
195                        f_result_.resize(n_order);
196                        //
197                        for(size_t k = 0; k < n_order; k++)
198                        {       f_left_[k].resize(nr_left_, n_middle_);
199                                f_right_[k].resize(n_middle_, nc_right_);
200                                f_result_[k].resize(nr_left_, nc_right_);
201                        }
202                }
203                // -------------------------------------------------------------------
204                // unpack tx into f_left and f_right
205                for(size_t k = 0; k < n_order; k++)
206                {       // unpack left values for this order
207                        for(size_t i = 0; i < n_left; i++)
208                                f_left_[k].data()[i] = tx[ i * n_order + k ];
209                        //
210                        // unpack right values for this order
211                        for(size_t i = 0; i < n_right; i++)
212                                f_right_[k].data()[i] = tx[ (i + n_left) * n_order + k ];
213                }
214                // -------------------------------------------------------------------
215                // result for each order
216                for(size_t k = 0; k < n_order; k++)
217                {       // result[k] = sum_ell left[ell] * right[k-ell]
218                        f_result_[k] = matrix::Zero(nr_left_, nc_right_);
219                        for(size_t ell = 0; ell <= k; ell++)
220                                f_result_[k] += f_left_[ell] * f_right_[k-ell];
221                }
222                // -------------------------------------------------------------------
223                // pack result_ into ty
224                for(size_t k = 0; k < n_order; k++)
225                {       for(size_t i = 0; i < n_result; i++)
226                                ty[ i * n_order + k ] = f_result_[k].data()[i];
227                }
228
229                // check if we are compute vy
230                if( vx.size() == 0 )
231                        return true;
232                // ------------------------------------------------------------------
233                // compute variable information for y; i.e., vy
234                // (note that the constant zero times a variable is a constant)
235                scalar zero(0.0);
236                assert( n_order == 1 );
237                for(size_t i = 0; i < nr_left_; i++)
238                {       for(size_t j = 0; j < nc_right_; j++)
239                        {       bool var = false;
240                                for(size_t ell = 0; ell < n_middle_; ell++)
241                                {       size_t index   = i * n_middle_ + ell;
242                                        bool var_left  = vx[index];
243                                        bool nz_left   = var_left | (f_left_[0](i, ell) != zero);
244                                        index          = nr_left_ * n_middle_;
245                                        index         += ell * nc_right_ + j;
246                                        bool var_right = vx[index];
247                                        bool nz_right  = var_right | (f_right_[0](ell, j) != zero);
248                                        var |= var_left & nz_right;
249                                        var |= nz_left  & var_right;
250                                }
251                                size_t index = i * nc_right_ + j;
252                                vy[index]    = var;
253                        }
254                }
255                return true;
256        }
257/* %$$
258$subhead reverse$$
259$srccode%cpp% */
260        // reverse mode routine called by CppAD
261        virtual bool reverse(
262                // highest order Taylor coefficient that we are computing deritive of
263                size_t                     q ,
264                // forward mode Taylor coefficients for x variables
265                const CppAD::vector<double>&     tx ,
266                // forward mode Taylor coefficients for y variables
267                const CppAD::vector<double>&     ty ,
268                // upon return, derivative of G[ F[ {x_j^k} ] ] w.r.t {x_j^k}
269                CppAD::vector<double>&           px ,
270                // derivative of G[ {y_i^k} ] w.r.t. {y_i^k}
271                const CppAD::vector<double>&     py
272        )
273        {       size_t n_order = q + 1;
274                assert( nx_ * n_order == tx.size() );
275                assert( ny_ * n_order == ty.size() );
276                assert( px.size() == tx.size() );
277                assert( py.size() == ty.size() );
278                //
279                size_t n_left   = nr_left_ * n_middle_;
280                size_t n_right  = n_middle_ * nc_right_;
281                size_t n_result = nr_left_ * nc_right_;
282                assert( n_left + n_right == nx_ );
283                assert( n_result == ny_ );
284                //
285                // -------------------------------------------------------------------
286                // make sure r_left_, r_right_, and r_result_ are large enough
287                assert( r_left_.size() == r_right_.size() );
288                assert( r_left_.size() == r_result_.size() );
289                if( r_left_.size() < n_order )
290                {       r_left_.resize(n_order);
291                        r_right_.resize(n_order);
292                        r_result_.resize(n_order);
293                        //
294                        for(size_t k = 0; k < n_order; k++)
295                        {       r_left_[k].resize(nr_left_, n_middle_);
296                                r_right_[k].resize(n_middle_, nc_right_);
297                                r_result_[k].resize(nr_left_, nc_right_);
298                        }
299                }
300                // -------------------------------------------------------------------
301                // unpack tx into f_left and f_right
302                for(size_t k = 0; k < n_order; k++)
303                {       // unpack left values for this order
304                        for(size_t i = 0; i < n_left; i++)
305                                f_left_[k].data()[i] = tx[ i * n_order + k ];
306                        //
307                        // unpack right values for this order
308                        for(size_t i = 0; i < n_right; i++)
309                                f_right_[k].data()[i] = tx[ (i + n_left) * n_order + k ];
310                }
311                // -------------------------------------------------------------------
312                // unpack result_ from py
313                for(size_t k = 0; k < n_order; k++)
314                {       for(size_t i = 0; i < n_result; i++)
315                                r_result_[k].data()[i] = py[ i * n_order + k ];
316                }
317                // -------------------------------------------------------------------
318                // initialize r_left_ and r_right_ as zero
319                for(size_t k = 0; k < n_order; k++)
320                {       r_left_[k]   = matrix::Zero(nr_left_, n_middle_);
321                        r_right_[k]  = matrix::Zero(n_middle_, nc_right_);
322                }
323                // -------------------------------------------------------------------
324                // matrix reverse mode calculation
325                for(size_t k1 = n_order; k1 > 0; k1--)
326                {       size_t k = k1 - 1;
327                        for(size_t ell = 0; ell <= k; ell++)
328                        {       // nr x nm       = nr x nc      * nc * nm
329                                r_left_[ell]    += r_result_[k] * f_right_[k-ell].transpose();
330                                // nm x nc       = nm x nr * nr * nc
331                                r_right_[k-ell] += f_left_[ell].transpose() * r_result_[k];
332                        }
333                }
334                // -------------------------------------------------------------------
335                // pack r_left and r_right int px
336                for(size_t k = 0; k < n_order; k++)
337                {       // pack left values for this order
338                        for(size_t i = 0; i < n_left; i++)
339                                px[ i * n_order + k ] = r_left_[k].data()[i];
340                        //
341                        // pack right values for this order
342                        for(size_t i = 0; i < n_right; i++)
343                                px[ (i + n_left) * n_order + k] = r_right_[k].data()[i];
344                }
345                //
346                return true;
347        }
348/* %$$
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/* %$$
528$head End Class Definition$$
529$srccode%cpp% */
530}; // End of atomic_eigen_mat_mul class
531
532}  // END_EMPTY_NAMESPACE
533/* %$$
534$$ $comment end nospell$$
535$end
536*/
537
538
539# endif
Note: See TracBrowser for help on using the repository browser.