source: trunk/example/sparse/rc_sparsity.cpp @ 3900

Last change on this file since 3900 was 3900, checked in by bradbell, 3 years ago

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

commit 3c93ecbd204908bd3b5e9cdb5c90824218a49356
Author: Brad Bell <bradbell@…>
Date: Fri Mar 10 06:00:33 2017 -0700

Editis replated to deprectated features removed on 2017-02-10.
run_cmake.sh: remove option to include these featurues.
configure.ac: group --with-deprecated with other flags that have been removed.

commit 1b68b229aba41c839582b678ac53d891e7893daa
Author: Brad Bell <bradbell@…>
Date: Fri Mar 10 05:16:38 2017 -0700

  1. Add sizing constructor to sparse_rc.
  2. Advanced to cppad-20170310.
File size: 9.4 KB
Line 
1/* --------------------------------------------------------------------------
2CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
3
4CppAD is distributed under multiple licenses. This distribution is under
5the terms of the
6                    Eclipse Public License Version 1.0.
7
8A copy of this license is included in the COPYING file of this distribution.
9Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
10-------------------------------------------------------------------------- */
11
12/*
13$begin rc_sparsity.cpp$$
14$spell
15        Bool
16        Jacobians
17$$
18
19$section Preferred Sparsity Patterns: Row and Column Indices: Example and Test$$
20
21$head Purpose$$
22This example show how to use row and column index sparsity patterns
23$cref sparse_rc$$ to compute sparse Jacobians and Hessians.
24This became the preferred way to represent sparsity on
25$cref/2017-02-09/whats_new_17/02-09/$$.
26
27$code
28$srcfile%example/sparse/rc_sparsity.cpp%0%// BEGIN C++%// END C++%1%$$
29$$
30
31$end
32*/
33// BEGIN C++
34# include <cppad/cppad.hpp>
35namespace {
36        using CppAD::sparse_rc;
37        using CppAD::sparse_rcv;
38        using CppAD::NearEqual;
39        //
40        typedef CPPAD_TESTVECTOR(bool)                b_vector;
41        typedef CPPAD_TESTVECTOR(size_t)              s_vector;
42        typedef CPPAD_TESTVECTOR(double)              d_vector;
43        typedef CPPAD_TESTVECTOR( CppAD::AD<double> ) a_vector;
44        //
45        double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
46        // -----------------------------------------------------------------------
47        // function f(x) that we are computing sparse results for
48        // -----------------------------------------------------------------------
49        a_vector fun(const a_vector& x)
50        {       size_t n  = x.size();
51                a_vector ret(n + 1);
52                for(size_t i = 0; i < n; i++)
53                {       size_t j = (i + 1) % n;
54                        ret[i]     = x[i] * x[i] * x[j];
55                }
56                ret[n] = 0.0;
57                return ret;
58        }
59        // -----------------------------------------------------------------------
60        // Jacobian
61        // -----------------------------------------------------------------------
62        bool check_jac(
63                const d_vector&                       x      ,
64                const sparse_rcv<s_vector, d_vector>& subset )
65        {       bool ok  = true;
66                size_t n = x.size();
67                //
68                ok &= subset.nnz() == 2 * n;
69                const s_vector& row( subset.row() );
70                const s_vector& col( subset.col() );
71                const d_vector& val( subset.val() );
72                s_vector row_major = subset.row_major();
73                for(size_t i = 0; i < n; i++)
74                {       size_t j = (i + 1) % n;
75                        size_t k = 2 * i;
76                        //
77                        ok &= row[ row_major[k] ]   == i;
78                        ok &= row[ row_major[k+1] ] == i;
79                        //
80                        size_t ck  = col[ row_major[k] ];
81                        size_t ckp = col[ row_major[k+1] ];
82                        double vk  = val[ row_major[k] ];
83                        double vkp = val[ row_major[k+1] ];
84                        //
85                        // put diagonal element first
86                        if( j < i )
87                        {       std::swap(ck, ckp);
88                                std::swap(vk, vkp);
89                        }
90                        // diagonal element
91                        ok &= ck == i;
92                        ok &= NearEqual( vk, 2.0 * x[i] * x[j], eps99, eps99 );
93                        // off diagonal element
94                        ok &= ckp == j;
95                        ok &= NearEqual( vkp, x[i] * x[i], eps99, eps99 );
96                }
97                return ok;
98        }
99        // Use forward mode for Jacobian and sparsity pattern
100        bool forward_jac(CppAD::ADFun<double>& f)
101        {       bool ok = true;
102                size_t n = f.Domain();
103                //
104                // sparsity pattern for identity matrix
105                sparse_rc<s_vector> pattern_in(n, n, n);
106                for(size_t k = 0; k < n; k++)
107                        pattern_in.set(k, k, k);
108                //
109                // sparsity pattern for Jacobian
110                bool transpose     = false;
111                bool dependency    = false;
112                bool internal_bool = false;
113                sparse_rc<s_vector> pattern_out;
114                f.for_jac_sparsity(
115                        pattern_in, transpose, dependency, internal_bool, pattern_out
116                );
117                //
118                // compute entire Jacobian
119                size_t                         group_max = 1;
120                std::string                    coloring  = "cppad";
121                sparse_rcv<s_vector, d_vector> subset( pattern_out );
122                CppAD::sparse_jac_work         work;
123                d_vector x(n);
124                for(size_t j = 0; j < n; j++)
125                        x[j] = double(j + 2);
126                size_t n_sweep = f.sparse_jac_for(
127                        group_max, x, subset, pattern_out, coloring, work
128                );
129                //
130                // check Jacobian
131                ok &= check_jac(x, subset);
132                ok &= n_sweep == 2;
133                //
134                return ok;
135        }
136        // Use reverse mode for Jacobian and sparsity pattern
137        bool reverse_jac(CppAD::ADFun<double>& f)
138        {       bool ok = true;
139                size_t n = f.Domain();
140                size_t m = f.Range();
141                //
142                // sparsity pattern for identity matrix
143                sparse_rc<s_vector> pattern_in(m, m, m);
144                for(size_t k = 0; k < m; k++)
145                        pattern_in.set(k, k, k);
146                //
147                // sparsity pattern for Jacobian
148                bool transpose     = false;
149                bool dependency    = false;
150                bool internal_bool = false;
151                sparse_rc<s_vector> pattern_out;
152                f.rev_jac_sparsity(
153                        pattern_in, transpose, dependency, internal_bool, pattern_out
154                );
155                //
156                // compute entire Jacobian
157                std::string                    coloring  = "cppad";
158                sparse_rcv<s_vector, d_vector> subset( pattern_out );
159                CppAD::sparse_jac_work         work;
160                d_vector x(n);
161                for(size_t j = 0; j < n; j++)
162                        x[j] = double(j + 2);
163                size_t n_sweep = f.sparse_jac_rev(
164                        x, subset, pattern_out, coloring, work
165                );
166                //
167                // check Jacobian
168                ok &= check_jac(x, subset);
169                ok &= n_sweep == 2;
170                //
171                return ok;
172        }
173        // ------------------------------------------------------------------------
174        // Hessian
175        // ------------------------------------------------------------------------
176        bool check_hes(
177                size_t                                i      ,
178                const d_vector&                       x      ,
179                const sparse_rcv<s_vector, d_vector>& subset )
180        {       bool ok  = true;
181                size_t n = x.size();
182                size_t j = (i + 1) % n;
183                //
184                ok &= subset.nnz() == 3;
185                const s_vector& row( subset.row() );
186                const s_vector& col( subset.col() );
187                const d_vector& val( subset.val() );
188                s_vector row_major = subset.row_major();
189                //
190                double v0 = val[ row_major[0] ];
191                double v1 = val[ row_major[1] ];
192                double v2 = val[ row_major[2] ];
193                if( j < i )
194                {       ok &= row[ row_major[0] ] == j;
195                        ok &= col[ row_major[0] ] == i;
196                        ok &= NearEqual( v0, 2.0 * x[i], eps99, eps99 );
197                        //
198                        ok &= row[ row_major[1] ] == i;
199                        ok &= col[ row_major[1] ] == j;
200                        ok &= NearEqual( v1, 2.0 * x[i], eps99, eps99 );
201                        //
202                        ok &= row[ row_major[2] ] == i;
203                        ok &= col[ row_major[2] ] == i;
204                        ok &= NearEqual( v2, 2.0 * x[j], eps99, eps99 );
205                }
206                else
207                {       ok &= row[ row_major[0] ] == i;
208                        ok &= col[ row_major[0] ] == i;
209                        ok &= NearEqual( v0, 2.0 * x[j], eps99, eps99 );
210                        //
211                        ok &= row[ row_major[1] ] == i;
212                        ok &= col[ row_major[1] ] == j;
213                        ok &= NearEqual( v1, 2.0 * x[i], eps99, eps99 );
214                        //
215                        ok &= row[ row_major[2] ] == j;
216                        ok &= col[ row_major[2] ] == i;
217                        ok &= NearEqual( v2, 2.0 * x[i], eps99, eps99 );
218                }
219                return ok;
220        }
221        // Use forward mode for Hessian and sparsity pattern
222        bool forward_hes(CppAD::ADFun<double>& f)
223        {       bool ok = true;
224                size_t n = f.Domain();
225                size_t m = f.Range();
226                //
227                b_vector select_domain(n);
228                for(size_t j = 0; j < n; j++)
229                        select_domain[j] = true;
230                sparse_rc<s_vector> pattern_out;
231                //
232                for(size_t i = 0; i < m; i++)
233                {       // select i-th component of range
234                        b_vector select_range(m);
235                        d_vector w(m);
236                        for(size_t k = 0; k < m; k++)
237                        {       select_range[k] = k == i;
238                                w[k] = 0.0;
239                                if( k == i )
240                                        w[k] = 1.0;
241                        }
242                        //
243                        bool internal_bool = false;
244                        f.for_hes_sparsity(
245                                select_domain, select_range, internal_bool, pattern_out
246                        );
247                        //
248                        // compute Hessian for i-th component function
249                        std::string                    coloring  = "cppad.symmetric";
250                        sparse_rcv<s_vector, d_vector> subset( pattern_out );
251                        CppAD::sparse_hes_work         work;
252                        d_vector x(n);
253                        for(size_t j = 0; j < n; j++)
254                                x[j] = double(j + 2);
255                        size_t n_sweep = f.sparse_hes(
256                                x, w, subset, pattern_out, coloring, work
257                        );
258                        //
259                        // check Hessian
260                        if( i == n )
261                                ok &= subset.nnz() == 0;
262                        else
263                        {       ok &= check_hes(i, x, subset);
264                                ok &= n_sweep == 2;
265                        }
266                }
267                return ok;
268        }
269        // Use reverse mode for Hessian and sparsity pattern
270        bool reverse_hes(CppAD::ADFun<double>& f)
271        {       bool ok = true;
272                size_t n = f.Domain();
273                size_t m = f.Range();
274                //
275                // n by n identity matrix
276                sparse_rc<s_vector> pattern_in(n, n, n);
277                for(size_t j = 0; j < n; j++)
278                        pattern_in.set(j, j, j);
279                //
280                bool transpose     = false;
281                bool dependency    = false;
282                bool internal_bool = true;
283                sparse_rc<s_vector> pattern_out;
284                //
285                f.for_jac_sparsity(
286                        pattern_in, transpose, dependency, internal_bool, pattern_out
287                );
288                //
289                for(size_t i = 0; i < m; i++)
290                {       // select i-th component of range
291                        b_vector select_range(m);
292                        d_vector w(m);
293                        for(size_t k = 0; k < m; k++)
294                        {       select_range[k] = k == i;
295                                w[k] = 0.0;
296                                if( k == i )
297                                        w[k] = 1.0;
298                        }
299                        //
300                        f.rev_hes_sparsity(
301                                select_range, transpose, internal_bool, pattern_out
302                        );
303                        //
304                        // compute Hessian for i-th component function
305                        std::string                    coloring  = "cppad.symmetric";
306                        sparse_rcv<s_vector, d_vector> subset( pattern_out );
307                        CppAD::sparse_hes_work         work;
308                        d_vector x(n);
309                        for(size_t j = 0; j < n; j++)
310                                x[j] = double(j + 2);
311                        size_t n_sweep = f.sparse_hes(
312                                x, w, subset, pattern_out, coloring, work
313                        );
314                        //
315                        // check Hessian
316                        if( i == n )
317                                ok &= subset.nnz() == 0;
318                        else
319                        {       ok &= check_hes(i, x, subset);
320                                ok &= n_sweep == 2;
321                        }
322                }
323                return ok;
324        }
325}
326// driver for all of the cases above
327bool rc_sparsity(void)
328{       bool ok = true;
329        //
330        // record the funcion
331        size_t n = 20;
332        size_t m = n + 1;
333        a_vector x(n), y(m);
334        for(size_t j = 0; j < n; j++)
335                x[j] = CppAD::AD<double>(j+1);
336        CppAD::Independent(x);
337        y = fun(x);
338        CppAD::ADFun<double> f(x, y);
339        //
340        // run the example / tests
341        ok &= forward_jac(f);
342        ok &= reverse_jac(f);
343        ok &= forward_hes(f);
344        ok &= reverse_hes(f);
345        //
346        return ok;
347}
348// END C++
Note: See TracBrowser for help on using the repository browser.