1 # ifndef CPPAD_CORE_SUBGRAPH_JAC_REV_HPP
2 # define CPPAD_CORE_SUBGRAPH_JAC_REV_HPP
3 /* --------------------------------------------------------------------------
4 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-18 Bradley M. Bell
5
6 CppAD is distributed under the terms of the
7 Eclipse Public License Version 2.0.
8
9 This Source Code may also be made available under the following
10 Secondary License when the conditions for such availability set forth
11 in the Eclipse Public License, Version 2.0 are satisfied:
12 GNU General Public License, Version 2.0 or later.
13 ---------------------------------------------------------------------------- */
14 /*
15 $begin subgraph_jac_rev$$
16 $spell
17 Jacobians
18 Jacobian
19 Subgraphs
20 subgraph
21 jac
22 rcv
23 Taylor
24 rev
25 nr
26 nc
27 const
28 Bool
29 nnz
30 $$
31
32 $section Compute Sparse Jacobians Using Subgraphs$$
33
34 $head Syntax$$
35 $icode%f%.subgraph_jac_rev(%x%, %subset%)
36 %$$
37 $icode%f%.subgraph_jac_rev(
38 %select_domain%, %select_range%, %x%, %matrix_out%
39 )%$$
40
41 $head See Also$$
42 $cref/clear_subgraph/subgraph_reverse/clear_subgraph/$$.
43
44 $head Purpose$$
45 We use $latex F : \B{R}^n \rightarrow \B{R}^m$$ to denote the
46 function corresponding to $icode f$$.
47 Here $icode n$$ is the $cref/domain/seq_property/Domain/$$ size,
48 and $icode m$$ is the $cref/range/seq_property/Range/$$ size, or $icode f$$.
49 The syntax above takes advantage of sparsity when computing the Jacobian
50 $latex \[
51 J(x) = F^{(1)} (x)
52 \] $$
53 The first syntax requires one to know what which elements of the Jacobian
54 they want to compute.
55 The second syntax computes the sparsity pattern and the value
56 of the Jacobian at the same time.
57 If one only wants the sparsity pattern,
58 it should be faster to use $cref subgraph_sparsity$$.
59
60 $head Method$$
61 This routine uses a subgraph technique. To be specific,
62 for each dependent variable,
63 it creates a subgraph of the operation sequence
64 containing the variables that affect the dependent variable.
65 This avoids the overhead of performing set operations
66 that is inherent in other methods for computing sparsity patterns.
67
68 $head BaseVector$$
69 The type $icode BaseVector$$ is a $cref SimpleVector$$ class with
70 $cref/elements of type/SimpleVector/Elements of Specified Type/$$
71 $icode Base$$.
72
73 $head SizeVector$$
74 The type $icode SizeVector$$ is a $cref SimpleVector$$ class with
75 $cref/elements of type/SimpleVector/Elements of Specified Type/$$
76 $code size_t$$.
77
78 $head BoolVector$$
79 The type $icode BoolVector$$ is a $cref SimpleVector$$ class with
80 $cref/elements of type/SimpleVector/Elements of Specified Type/$$
81 $code bool$$.
82
83 $head f$$
84 This object has prototype
85 $codei%
86 ADFun<%Base%> %f%
87 %$$
88 Note that the Taylor coefficients stored in $icode f$$ are affected
89 by this operation; see
90 $cref/uses forward/sparse_jac/Uses Forward/$$ below.
91
92 $head x$$
93 This argument has prototype
94 $codei%
95 const %BaseVector%& %x%
96 %$$
97 It is the value of $icode x$$ at which we are computing the Jacobian.
98
99 $head Uses Forward$$
100 After each call to $cref Forward$$,
101 the object $icode f$$ contains the corresponding
102 $cref/Taylor coefficients/glossary/Taylor Coefficient/$$.
103 After a call to $code sparse_jac_forward$$ or $code sparse_jac_rev$$,
104 the zero order coefficients correspond to
105 $codei%
106 %f%.Forward(0, %x%)
107 %$$
108 All the other forward mode coefficients are unspecified.
109
110 $head subset$$
111 This argument has prototype
112 $codei%
113 sparse_rcv<%SizeVector%, %BaseVector%>& %subset%
114 %$$
115 Its row size is $icode%subset%.nr() == %m%$$,
116 and its column size is $icode%subset%.nc() == %n%$$.
117 It specifies which elements of the Jacobian are computed.
118 The input elements in its value vector
119 $icode%subset%.val()%$$ do not matter.
120 Upon return it contains the value of the corresponding elements
121 of the Jacobian.
122
123 $head select_domain$$
124 The argument $icode select_domain$$ has prototype
125 $codei%
126 const %BoolVector%& %select_domain%
127 %$$
128 It has size $latex n$$ and specifies which independent variables
129 to include.
130
131 $head select_range$$
132 The argument $icode select_range$$ has prototype
133 $codei%
134 const %BoolVector%& %select_range%
135 %$$
136 It has size $latex m$$ and specifies which components of the range
137 to include in the calculation.
138 A subgraph is built for each dependent variable and the selected set
139 of independent variables.
140
141 $head matrix_out$$
142 This argument has prototype
143 $codei%
144 sparse_rcv<%SizeVector%, %BaseVector%>& %matrix_out%
145 %$$
146 This input value of $icode matrix_out$$ does not matter.
147 Upon return $icode matrix_out$$ is
148 $cref/sparse matrix/sparse_rcv/$$ representation of $latex F^{(1)} (x)$$.
149 The matrix has $latex m$$ rows, $latex n$$ columns.
150 If $icode%select_domain%[%j%]%$$ is true,
151 $icode%select_range%[%i%]%$$ is true, and
152 $latex F_i (x)$$ depends on $latex x_j$$,
153 then the pair $latex (i, j)$$ is in $icode matrix_out$$.
154 For each $icode%k% = 0 , %...%, %matrix_out%.nnz()%$$, let
155 $codei%
156 %i% = %matrix_out%.row()[%k%]
157 %j% = %matrix_out%.col()[%k%]
158 %v% = %matrix_out%.val()[%k%]
159 %$$
160 It follows that the partial of $latex F_i (x)$$ with respect to
161 $latex x_j$$ is equal to $latex v$$.
162
163
164 $head Example$$
165 $children%
166 example/sparse/subgraph_jac_rev.cpp%
167 example/sparse/subgraph_hes2jac.cpp
168 %$$
169 The files $cref subgraph_jac_rev.cpp$$ and $cref subgraph_hes2jac.cpp$$
170 are examples and tests using $code subgraph_jac_rev$$.
171 They returns $code true$$ for success and $code false$$ for failure.
172
173 $end
174 -----------------------------------------------------------------------------
175 */
176 # include <cppad/core/ad_fun.hpp>
177 # include <cppad/local/subgraph/info.hpp>
178
179 namespace CppAD { // BEGIN_CPPAD_NAMESPACE
180
181 /*!
182 Subgraph sparsity patterns.
183
184 \tparam Base
185 is the base type for this recording.
186
187 \tparam SizeVector
188 is the simple vector with elements of type size_t that is used for
189 row, column index sparsity patterns.
190
191 \tparam BaseVector
192 a simple vector class with elements of type Base.
193
194 \param x
195 a vector of length n, the number of independent variables in f
196 (this ADFun object).
197
198 \param subset
199 specifices the subset of the sparsity pattern where the Jacobian is evaluated.
200 subset.nr() == m,
201 subset.nc() == n.
202 */
203 template <class Base, class RecBase>
204 template <class SizeVector, class BaseVector>
subgraph_jac_rev(const BaseVector & x,sparse_rcv<SizeVector,BaseVector> & subset)205 void ADFun<Base,RecBase>::subgraph_jac_rev(
206 const BaseVector& x ,
207 sparse_rcv<SizeVector, BaseVector>& subset )
208 { size_t m = Range();
209 size_t n = Domain();
210 //
211 CPPAD_ASSERT_KNOWN(
212 subset.nr() == m,
213 "subgraph_jac_rev: subset.nr() not equal range dimension for f"
214 );
215 CPPAD_ASSERT_KNOWN(
216 subset.nc() == n,
217 "subgraph_jac_rev: subset.nc() not equal domain dimension for f"
218 );
219 //
220 // point at which we are evaluating Jacobian
221 Forward(0, x);
222 //
223 // nnz and row, column, and row_major vectors for subset
224 size_t nnz = subset.nnz();
225 const SizeVector& row( subset.row() );
226 const SizeVector& col( subset.col() );
227 SizeVector row_major = subset.row_major();
228 //
229 // determine set of independent variabels
230 local::pod_vector<bool> select_domain(n);
231 for(size_t j = 0; j < n; j++)
232 select_domain[j] = false;
233 for(size_t k = 0; k < nnz; k++)
234 select_domain[ col[k] ] = true;
235 //
236 // initialize reverse mode computation on subgraphs
237 subgraph_reverse(select_domain);
238 //
239 // memory used to hold subgraph_reverse results
240 BaseVector dw;
241 SizeVector dw_col;
242 //
243 // initialize index in row_major
244 size_t k = 0;
245 Base zero(0);
246 while(k < nnz )
247 { size_t q = 1;
248 size_t i_dep = row[ row_major[k] ];
249 size_t i_ind = col[ row_major[k] ];
250 size_t ell = i_dep;
251 subgraph_reverse(q, ell, dw_col, dw);
252 //
253 size_t c = 0;
254 while( i_dep == ell )
255 { // row numbers match
256 //
257 // advance c to possible match with column i_ind
258 while( c < size_t( dw_col.size() ) && dw_col[c] < i_ind )
259 ++c;
260 //
261 // check for match with i_ind
262 if( i_ind == dw_col[c] )
263 subset.set( row_major[k], dw[i_ind] );
264 else
265 subset.set( row_major[k], zero);
266 //
267 // advance to next (i_dep, i_ind)
268 ++k;
269 if( k == nnz )
270 { i_dep = m;
271 i_ind = n;
272 }
273 else
274 { i_dep = row[ row_major[k] ];
275 i_ind = col[ row_major[k] ];
276 }
277 }
278 }
279 return;
280 }
281 template <class Base, class RecBase>
282 template <class BoolVector, class SizeVector, class BaseVector>
subgraph_jac_rev(const BoolVector & select_domain,const BoolVector & select_range,const BaseVector & x,sparse_rcv<SizeVector,BaseVector> & matrix_out)283 void ADFun<Base,RecBase>::subgraph_jac_rev(
284 const BoolVector& select_domain ,
285 const BoolVector& select_range ,
286 const BaseVector& x ,
287 sparse_rcv<SizeVector, BaseVector>& matrix_out )
288 { size_t m = Range();
289 size_t n = Domain();
290 //
291 // point at which we are evaluating Jacobian
292 Forward(0, x);
293 //
294 // nnz and row, column, and row_major vectors for subset
295 local::pod_vector<size_t> row_out;
296 local::pod_vector<size_t> col_out;
297 local::pod_vector_maybe<Base> val_out;
298 //
299 // initialize reverse mode computation on subgraphs
300 subgraph_reverse(select_domain);
301 //
302 // memory used to hold subgraph_reverse results
303 BaseVector dw;
304 SizeVector col;
305 //
306 // loop through selected independent variables
307 for(size_t i = 0; i < m; ++i) if( select_range[i] )
308 { // compute Jacobian and sparsity for this dependent variable
309 size_t q = 1;
310 subgraph_reverse(q, i, col, dw);
311 CPPAD_ASSERT_UNKNOWN( size_t( dw.size() ) == n );
312 //
313 // offset for this dependent variable
314 size_t index = row_out.size();
315 CPPAD_ASSERT_UNKNOWN( col_out.size() == index );
316 CPPAD_ASSERT_UNKNOWN( val_out.size() == index );
317 //
318 // extend vectors to hold results for this dependent variable
319 size_t col_size = size_t( col.size() );
320 row_out.extend( col_size );
321 col_out.extend( col_size );
322 val_out.extend( col_size );
323 //
324 // store results for this dependent variable
325 for(size_t c = 0; c < col_size; ++c)
326 { row_out[index + c] = i;
327 col_out[index + c] = col[c];
328 val_out[index + c] = dw[ col[c] ];
329 }
330 }
331 //
332 // create sparsity pattern corresponding to row_out, col_out
333 size_t nr = m;
334 size_t nc = n;
335 size_t nnz = row_out.size();
336 sparse_rc<SizeVector> pattern(nr, nc, nnz);
337 for(size_t k = 0; k < nnz; ++k)
338 pattern.set(k, row_out[k], col_out[k]);
339 //
340 // create sparse matrix
341 sparse_rcv<SizeVector, BaseVector> matrix(pattern);
342 for(size_t k = 0; k < nnz; ++k)
343 matrix.set(k, val_out[k]);
344 //
345 // return matrix
346 matrix_out = matrix;
347 //
348 return;
349 }
350 } // END_CPPAD_NAMESPACE
351 # endif
352