source: trunk/example/atomic/eigen_mat_mul.cpp @ 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: 8.3 KB
Line 
1// $Id$
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-16 Bradley M. Bell
4
5CppAD is distributed under multiple licenses. This distribution is under
6the terms of the
7                    Eclipse Public License Version 1.0.
8
9A copy of this license is included in the COPYING file of this distribution.
10Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
11-------------------------------------------------------------------------- */
12
13/*
14$begin atomic_eigen_mat_mul.cpp$$
15$spell
16        mul
17        Eigen
18$$
19
20$section  Atomic Eigen Matrix Multiply: Example and Test$$
21
22$head Under Construction$$
23This example is under construction. So far only forward mode has been
24implemented.
25
26$head Description$$
27The $cref ADFun$$ function object $icode f$$ for this example is
28$latex \[
29f(x) =
30\left( \begin{array}{cc}
31        0   & 0 \\
32        1   & 2 \\
33        x_0 & x_1
34\end{array} \right)
35\left( \begin{array}{c}
36        x_0 \\
37        x_1
38\end{array} \right)
39=
40\left( \begin{array}{c}
41        0 \\
42        x_0 + 2 x_1 \\
43        x_0 x_0 + x_1 x_1 )
44\end{array} \right)
45\] $$
46
47$children%
48        cppad/example/eigen_mat_mul.hpp
49%$$
50
51$head Class Definition$$
52This example uses the file $cref atomic_eigen_mat_mul.hpp$$
53which defines matrix multiply as a $cref atomic_base$$ operation.
54
55$nospell
56
57$head Use Atomic Function$$
58$srccode%cpp% */
59# include <cppad/cppad.hpp>
60# include <cppad/example/eigen_mat_mul.hpp>
61
62namespace {
63        typedef double            scalar;
64        typedef CppAD::AD<scalar> ad_scalar;
65        typedef typename atomic_eigen_mat_mul<scalar>::matrix     matrix;
66        typedef typename atomic_eigen_mat_mul<scalar>::ad_matrix  ad_matrix;
67
68        // use atomic operation to multiply two AD matrices
69        ad_matrix matrix_multiply(
70                atomic_eigen_mat_mul<scalar>& mat_mul ,
71                const ad_matrix&              left    ,
72                const ad_matrix&              right   )
73        {       size_t nr_left   = size_t( left.rows() );
74                size_t n_middle    = size_t( left.cols() );
75                size_t nc_right  = size_t ( right.cols() );
76                assert( size_t( right.rows() ) == n_middle );
77
78                // packed version of left and right
79                size_t nx = (nr_left + nc_right) * n_middle;
80                CPPAD_TESTVECTOR(ad_scalar) packed_arg(nx);
81                mat_mul.pack(packed_arg, left, right);
82
83                // packed version of result = left * right
84                size_t ny = nr_left * nc_right;
85                CPPAD_TESTVECTOR(ad_scalar) packed_result(ny);
86                mat_mul(packed_arg, packed_result);
87
88                // result matrix
89                ad_matrix result(nr_left, nc_right);
90                mat_mul.unpack(packed_result, result);
91
92                return result;
93        }
94
95}
96
97bool eigen_mat_mul(void)
98{       bool ok    = true;
99        scalar eps = 10. * std::numeric_limits<scalar>::epsilon();
100        using CppAD::NearEqual;
101        //
102/* %$$
103$subhead Constructor$$
104$srccode%cpp% */
105        // -------------------------------------------------------------------
106        // object that multiplies a 3x2 matrix times a 2x1 matrix
107        size_t nr_left  = 3;
108        size_t n_middle   = 2;
109        size_t nc_right = 1;
110        atomic_eigen_mat_mul<scalar> mat_mul(nr_left, n_middle, nc_right);
111        // -------------------------------------------------------------------
112        // declare independent variable vector x
113        size_t n = 2;
114        CPPAD_TESTVECTOR(ad_scalar) ad_x(n);
115        for(size_t j = 0; j < n; j++)
116                ad_x[j] = ad_scalar(j);
117        CppAD::Independent(ad_x);
118        // -------------------------------------------------------------------
119        //        [ 0     0    ]
120        // left = [ 1     2    ]
121        //        [ x[0]  x[1] ]
122        ad_matrix ad_left(nr_left, n_middle);
123        ad_left(0, 0) = ad_scalar(0.0);
124        ad_left(0, 1) = ad_scalar(0.0);
125        ad_left(1, 0) = ad_scalar(1.0);
126        ad_left(1, 1) = ad_scalar(2.0);
127        ad_left(2, 0) = ad_x[0];
128        ad_left(2, 1) = ad_x[1];
129        // -------------------------------------------------------------------
130        // right = [ x[0] , x[1] ]^T
131        ad_matrix ad_right(n_middle, nc_right);
132        ad_right(0, 0) = ad_x[0];
133        ad_right(1, 0) = ad_x[1];
134        // -------------------------------------------------------------------
135        // use atomic operation to multiply left * right
136        ad_matrix ad_result = matrix_multiply(mat_mul, ad_left, ad_right);
137        // -------------------------------------------------------------------
138        // check that first component of result is a parameter
139        // and the other components are varaibles.
140        ok &= Parameter( ad_result(0, 0) );
141        ok &= Variable(  ad_result(1, 0) );
142        ok &= Variable(  ad_result(2, 0) );
143        // -------------------------------------------------------------------
144        // declare the dependent variable vector y
145        size_t m = 3;
146        CPPAD_TESTVECTOR(ad_scalar) ad_y(m);
147        for(size_t i = 0; i < m; i++)
148                ad_y[i] = ad_result(i, 0);
149        CppAD::ADFun<scalar> f(ad_x, ad_y);
150        // -------------------------------------------------------------------
151        // check zero order forward mode
152        CPPAD_TESTVECTOR(scalar) x(n), y(m);
153        for(size_t i = 0; i < n; i++)
154                x[i] = scalar(i + 2);
155        y   = f.Forward(0, x);
156        ok &= NearEqual(y[0], 0.0,                       eps, eps);
157        ok &= NearEqual(y[1], x[0] + 2.0 * x[1],         eps, eps);
158        ok &= NearEqual(y[2], x[0] * x[0] + x[1] * x[1], eps, eps);
159        // -------------------------------------------------------------------
160        // check first order forward mode
161        CPPAD_TESTVECTOR(scalar) x1(n), y1(m);
162        x1[0] = 1.0;
163        x1[1] = 0.0;
164        y1    = f.Forward(1, x1);
165        ok   &= NearEqual(y1[0], 0.0,        eps, eps);
166        ok   &= NearEqual(y1[1], 1.0,        eps, eps);
167        ok   &= NearEqual(y1[2], 2.0 * x[0], eps, eps);
168        x1[0] = 0.0;
169        x1[1] = 1.0;
170        y1    = f.Forward(1, x1);
171        ok   &= NearEqual(y1[0], 0.0,        eps, eps);
172        ok   &= NearEqual(y1[1], 2.0,        eps, eps);
173        ok   &= NearEqual(y1[2], 2.0 * x[1], eps, eps);
174        // -------------------------------------------------------------------
175        // check second order forward mode
176        CPPAD_TESTVECTOR(scalar) x2(n), y2(m);
177        x2[0] = 0.0;
178        x2[1] = 0.0;
179        y2    = f.Forward(2, x2);
180        ok   &= NearEqual(y2[0], 0.0, eps, eps);
181        ok   &= NearEqual(y2[1], 0.0, eps, eps);
182        ok   &= NearEqual(y2[2], 1.0, eps, eps); // 1/2 * f_1''(x)
183        // -------------------------------------------------------------------
184        // check first order reverse mode
185        CPPAD_TESTVECTOR(scalar) w(m), d1w(n);
186        w[0]  = 0.0;
187        w[1]  = 1.0;
188        w[2]  = 0.0;
189        d1w   = f.Reverse(1, w);
190        ok   &= NearEqual(d1w[0], 1.0, eps, eps);
191        ok   &= NearEqual(d1w[1], 2.0, eps, eps);
192        w[0]  = 0.0;
193        w[1]  = 0.0;
194        w[2]  = 1.0;
195        d1w   = f.Reverse(1, w);
196        ok   &= NearEqual(d1w[0], 2.0 * x[0], eps, eps);
197        ok   &= NearEqual(d1w[1], 2.0 * x[1], eps, eps);
198        // -------------------------------------------------------------------
199        // check second order reverse mode
200        CPPAD_TESTVECTOR(scalar) d2w(2 * n);
201        d2w   = f.Reverse(2, w);
202        ok   &= NearEqual(d2w[0 * 2 + 0], 2.0 * x[0], eps, eps);
203        ok   &= NearEqual(d2w[1 * 2 + 0], 2.0 * x[1], eps, eps);
204        // partial f_1 w.r.t x_1, x_0
205        ok   &= NearEqual(d2w[0 * 2 + 1], 0.0,        eps, eps);
206        // partial f_1 w.r.t x_1, x_1
207        ok   &= NearEqual(d2w[1 * 2 + 1], 2.0,        eps, eps);
208        // -------------------------------------------------------------------
209        // check forward Jacobian sparsity
210        CPPAD_TESTVECTOR( std::set<size_t> ) r(n), s(m);
211        std::set<size_t> check_set;
212        for(size_t j = 0; j < n; j++)
213                r[j].insert(j);
214        s      = f.ForSparseJac(n, r);
215        check_set.clear();
216        ok    &= s[0] == check_set;
217        check_set.insert(0);
218        check_set.insert(1);
219        ok    &= s[1] == check_set;
220        ok    &= s[2] == check_set;
221        // -------------------------------------------------------------------
222        // check reverse Jacobian sparsity
223        r.resize(m);
224        for(size_t i = 0; i < m; i++)
225                r[i].insert(i);
226        s  = f.RevSparseJac(m, r);
227        check_set.clear();
228        ok    &= s[0] == check_set;
229        check_set.insert(0);
230        check_set.insert(1);
231        ok    &= s[1] == check_set;
232        ok    &= s[2] == check_set;
233        // -------------------------------------------------------------------
234        // check forward Hessian sparsity for f_2 (x)
235        CPPAD_TESTVECTOR( std::set<size_t> ) r2(1), s2(1), h(n);
236        for(size_t j = 0; j < n; j++)
237                r2[0].insert(j);
238        s2[0].clear();
239        s2[0].insert(2);
240        h = f.ForSparseHes(r2, s2);
241        check_set.clear();
242        check_set.insert(0);
243        ok &= h[0] == check_set;
244        check_set.clear();
245        check_set.insert(1);
246        ok &= h[1] == check_set;
247        // -------------------------------------------------------------------
248        // check reverse Hessian sparsity for f_2 (x)
249        CPPAD_TESTVECTOR( std::set<size_t> ) s3(1);
250        s3[0].clear();
251        s3[0].insert(2);
252        h = f.RevSparseHes(n, s3);
253        check_set.clear();
254        check_set.insert(0);
255        ok &= h[0] == check_set;
256        check_set.clear();
257        check_set.insert(1);
258        ok &= h[1] == check_set;
259        // -------------------------------------------------------------------
260        return ok;
261}
262/* %$$
263$$ $comment end nospell$$
264$end
265*/
Note: See TracBrowser for help on using the repository browser.