1 /* --------------------------------------------------------------------------
2 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-20 Bradley M. Bell
3
4 CppAD is distributed under the terms of the
5 Eclipse Public License Version 2.0.
6
7 This Source Code may also be made available under the following
8 Secondary License when the conditions for such availability set forth
9 in the Eclipse Public License, Version 2.0 are satisfied:
10 GNU General Public License, Version 2.0 or later.
11 ---------------------------------------------------------------------------- */
12 /*
13 $begin cppadcg_sparse_jacobian.cpp$$
14 $spell
15 cppadcg
16 jacobian
17 Cpp
18 $$
19
20 $section Cppadcg Speed: Sparse Jacobian$$
21
22 $head Specifications$$
23 See $cref link_sparse_jacobian$$.
24
25 $head PASS_SPARSE_JACOBIAN_TO_CODE_GEN$$
26 If this is one, the sparse Jacobian is the is the function passed
27 to CppADCodeGen, In this case, the $code code_gen_fun$$
28 $cref/function/code_gen_fun/Syntax/function/$$ is used to calculate
29 the sparse Jacobian.
30 Otherwise, this flag is zero and the original function passed
31 to CppADCodeGen. In this case, the $code code_gen_fun$$
32 $cref/sparse_jacobian/code_gen_fun/Syntax/sparse_jacobian/$$
33 is used to calculate the sparse Jacobian.
34 $srccode%cpp% */
35 # define PASS_SPARSE_JACOBIAN_TO_CODE_GEN 1
36 /* %$$
37
38
39 $head Implementation$$
40 $srccode%cpp% */
41 # include <cppad/speed/uniform_01.hpp>
42 # include <cppad/utility/vector.hpp>
43 # include <cppad/speed/sparse_jac_fun.hpp>
44 # include <cppad/example/code_gen_fun.hpp>
45
46 # include <map>
47 extern std::map<std::string, bool> global_option;
48
49 namespace {
50 // -----------------------------------------------------------------------
51 // typedefs
52 typedef CppAD::cg::CG<double> c_double;
53 typedef CppAD::AD<c_double> ac_double;
54 typedef CppAD::vector<bool> b_vector;
55 typedef CppAD::vector<size_t> s_vector;
56 typedef CppAD::vector<double> d_vector;
57 typedef CppAD::vector<ac_double> ac_vector;
58 typedef CppAD::sparse_rc<s_vector> sparsity;
59 // ------------------------------------------------------------------------
60 # if PASS_SPARSE_JACOBIAN_TO_CODE_GEN
61 // calc_sparsity
calc_sparsity(CppAD::sparse_rc<s_vector> & pattern,CppAD::ADFun<c_double> & f)62 void calc_sparsity(
63 CppAD::sparse_rc<s_vector>& pattern ,
64 CppAD::ADFun<c_double>& f )
65 { bool reverse = global_option["revsparsity"];
66 bool transpose = false;
67 bool internal_bool = global_option["boolsparsity"];
68 bool dependency = false;
69 bool subgraph = global_option["subsparsity"];
70 size_t n = f.Domain();
71 size_t m = f.Range();
72 if( subgraph )
73 { b_vector select_domain(n), select_range(m);
74 for(size_t j = 0; j < n; ++j)
75 select_domain[j] = true;
76 for(size_t i = 0; i < m; ++i)
77 select_range[i] = true;
78 f.subgraph_sparsity(
79 select_domain, select_range, transpose, pattern
80 );
81 }
82 else
83 { size_t q = n;
84 if( reverse )
85 q = m;
86 //
87 CppAD::sparse_rc<s_vector> identity;
88 identity.resize(q, q, q);
89 for(size_t k = 0; k < q; k++)
90 identity.set(k, k, k);
91 //
92 if( reverse )
93 { f.rev_jac_sparsity(
94 identity, transpose, dependency, internal_bool, pattern
95 );
96 }
97 else
98 { f.for_jac_sparsity(
99 identity, transpose, dependency, internal_bool, pattern
100 );
101 }
102 }
103 }
104 # endif // PASS_SPARSE_JACOBIAN_TO_CODE_GEN
105 // -------------------------------------------------------------------------
106 // setup
setup(size_t size,const s_vector & row,const s_vector & col,size_t & n_color,code_gen_fun & fun,s_vector & row_major)107 void setup(
108 // inputs
109 size_t size ,
110 const s_vector& row ,
111 const s_vector& col ,
112 // outputs
113 size_t& n_color ,
114 code_gen_fun& fun ,
115 s_vector& row_major )
116 { // optimization options
117 std::string optimize_options =
118 "no_conditional_skip no_compare_op no_print_for_op";
119 //
120 // independent variable vector
121 size_t nc = size;
122 ac_vector ac_x(nc);
123 //
124 // dependent variable vector
125 size_t nr = 2 * nc;
126 ac_vector ac_y(nr);
127 //
128 // values of independent variables do not matter
129 for(size_t j = 0; j < nc; j++)
130 ac_x[j] = ac_double( double(j) / double(nc) );
131 //
132 // declare independent variables
133 size_t abort_op_index = 0;
134 bool record_compare = false;
135 CppAD::Independent(ac_x, abort_op_index, record_compare);
136 //
137 // AD computation of f(x) (order zero derivative is function value)
138 size_t order = 0;
139 CppAD::sparse_jac_fun<ac_double>(nr, nc, ac_x, row, col, order, ac_y);
140 //
141 // create function object f : x -> y
142 CppAD::ADFun<c_double> c_f;
143 CppAD::ADFun<ac_double, c_double> ac_f;
144 c_f.Dependent(ac_x, ac_y);
145 if( global_option["optimize"] )
146 c_f.optimize(optimize_options);
147 //
148 // number of non-zeros in sparsity pattern for Jacobian
149 # if ! PASS_SPARSE_JACOBIAN_TO_CODE_GEN
150 // set fun
151 code_gen_fun::evaluation_enum eval_jac = code_gen_fun::sparse_enum;
152 code_gen_fun f_tmp("sparse_jacobian", c_f, eval_jac);
153 fun.swap(f_tmp);
154 //
155 // set row_major
156 d_vector x(nc);
157 CppAD::uniform_01(nc, x);
158 CppAD::sparse_rcv<s_vector, d_vector> Jrcv = fun.sparse_jacobian(x);
159 row_major = Jrcv.row_major();
160 # ifndef NDEBUG
161 size_t nnz = row.size();
162 CPPAD_ASSERT_UNKNOWN( row_major.size() == nnz );
163 for(size_t k = 0; k < nnz; ++k)
164 { size_t ell = row_major[k];
165 CPPAD_ASSERT_UNKNOWN(
166 Jrcv.row()[ell] == row[k] && Jrcv.col()[ell] == col[k]
167 );
168 }
169 # endif
170 //
171 # else // PASS_SPARSE_JACOBIAN_TO_CODE_GEN
172 //
173 // sparsity patttern for subset of Jacobian pattern that is evaluated
174 size_t nnz = row.size();
175 sparsity subset_pattern(nr, nc, nnz);
176 for(size_t k = 0; k < nnz; ++k)
177 subset_pattern.set(k, row[k], col[k]);
178 //
179 // spoarse matrix for subset of Jacobian that is evaluated
180 CppAD::sparse_rcv<s_vector, ac_vector> ac_subset( subset_pattern );
181 //
182 // coloring method
183 std::string coloring = "cppad";
184 # if CPPAD_HAS_COLPACK
185 if( global_option["colpack"] )
186 coloring = "colpack";
187 # endif
188 //
189 // maximum number of colors at once
190 size_t group_max = 1;
191 ac_f = c_f.base2ad();
192 //
193 // declare independent variables for jacobian computation
194 CppAD::Independent(ac_x, abort_op_index, record_compare);
195 //
196 if( global_option["subgraph"] )
197 { // use reverse mode because forward not yet implemented
198 ac_f.subgraph_jac_rev(ac_x, ac_subset);
199 n_color = 0;
200 }
201 else
202 { // calculate the Jacobian sparsity pattern for this function
203 sparsity pattern;
204 calc_sparsity(pattern, c_f);
205 //
206 // use forward mode to compute Jacobian
207 CppAD::sparse_jac_work work;
208 n_color = ac_f.sparse_jac_for(
209 group_max, ac_x, ac_subset, pattern, coloring, work
210 );
211 }
212 const ac_vector ac_val ( ac_subset.val() );
213 //
214 // create function g : x -> f'(x)
215 CppAD::ADFun<c_double> c_g;
216 c_g.Dependent(ac_x, ac_val);
217 if( global_option["optimize"] )
218 c_g.optimize(optimize_options);
219 code_gen_fun g_tmp("sparse_jacobian", c_g);
220 //
221 // set reture value
222 fun.swap(g_tmp);
223 # endif // PASS_SPARSE_JACOBIAN_TO_CODE_GEN
224 return;
225 }
226 }
227
link_sparse_jacobian(const std::string & job,size_t size,size_t repeat,size_t m,const CppAD::vector<size_t> & row,const CppAD::vector<size_t> & col,CppAD::vector<double> & x,CppAD::vector<double> & jacobian,size_t & n_color)228 bool link_sparse_jacobian(
229 const std::string& job ,
230 size_t size ,
231 size_t repeat ,
232 size_t m ,
233 const CppAD::vector<size_t>& row ,
234 const CppAD::vector<size_t>& col ,
235 CppAD::vector<double>& x ,
236 CppAD::vector<double>& jacobian ,
237 size_t& n_color )
238 { assert( x.size() == size );
239 assert( jacobian.size() == row.size() );
240 // --------------------------------------------------------------------
241 // check global options
242 const char* valid[] = {
243 "memory", "onetape", "optimize", "subgraph",
244 "boolsparsity", "revsparsity", "subsparsity"
245 # if CPPAD_HAS_COLPACK
246 , "colpack"
247 # endif
248 };
249 size_t n_valid = sizeof(valid) / sizeof(valid[0]);
250 typedef std::map<std::string, bool>::iterator iterator;
251 //
252 for(iterator itr=global_option.begin(); itr!=global_option.end(); ++itr)
253 { if( itr->second )
254 { bool ok = false;
255 for(size_t i = 0; i < n_valid; i++)
256 ok |= itr->first == valid[i];
257 if( ! ok )
258 return false;
259 }
260 }
261 if( global_option["subsparsity"] )
262 { if( global_option["boolsparisty"]
263 || global_option["revsparsity"]
264 || global_option["colpack"] )
265 return false;
266 }
267 // -----------------------------------------------------
268 // size corresponding to static_fun
269 static size_t static_size = 0;
270 //
271 // function object mapping x to f'(x)
272 static code_gen_fun static_fun;
273 //
274 // row_major order for Jrcv
275 static s_vector static_row_major;
276 //
277 # if ! PASS_SPARSE_JACOBIAN_TO_CODE_GEN
278 // code gen value for sparse jacobian
279 CppAD::sparse_rcv<s_vector, d_vector> Jrcv;
280 # endif
281 //
282 // number of independent variables
283 size_t nx = size;
284 //
285 bool onetape = global_option["onetape"];
286 //
287 // default return value
288 n_color = 0;
289 // -----------------------------------------------------
290 if( job == "setup" )
291 { if( onetape )
292 { // sets n_color when ontape is true
293 setup(size, row, col, n_color, static_fun, static_row_major);
294 static_size = size;
295 }
296 else
297 { static_size = 0;
298 }
299 return true;
300 }
301 if( job == "teardown" )
302 { code_gen_fun f_tmp;
303 static_fun.swap(f_tmp);
304 static_row_major.clear();
305 //
306 static_size = 0;
307 return true;
308 }
309 // -----------------------------------------------------
310 CPPAD_ASSERT_UNKNOWN( job == "run" )
311 if( onetape ) while(repeat--)
312 { // use if before assert to vaoid warning that static_size is not used
313 if( size != static_size )
314 { CPPAD_ASSERT_UNKNOWN( size == static_size );
315 }
316
317 // get next x
318 CppAD::uniform_01(nx, x);
319
320 // evaluate the jacobian
321 # if PASS_SPARSE_JACOBIAN_TO_CODE_GEN
322 jacobian = static_fun(x);
323 # else
324 Jrcv = static_fun.sparse_jacobian(x);
325 CPPAD_ASSERT_UNKNOWN( Jrcv.nnz() == jacobian.size() );
326 for(size_t k = 0; k < row.size(); ++k)
327 jacobian[k] = Jrcv.val()[ static_row_major[k] ];
328 # endif
329 }
330 else while(repeat--)
331 { // sets n_color when ontape is false
332 setup(size, row, col, n_color, static_fun, static_row_major);
333 static_size = size;
334
335 // get next x
336 CppAD::uniform_01(nx, x);
337
338 // evaluate the jacobian
339 # if PASS_SPARSE_JACOBIAN_TO_CODE_GEN
340 jacobian = static_fun(x);
341 # else
342 Jrcv = static_fun.sparse_jacobian(x);
343 CPPAD_ASSERT_UNKNOWN( Jrcv.nnz() == jacobian.size() );
344 for(size_t k = 0; k < row.size(); ++k)
345 jacobian[k] = Jrcv.val()[ static_row_major[k] ];
346 # endif
347 }
348 return true;
349 }
350 /* %$$
351 $end
352 */
353