CppAD: A C++ Algorithmic Differentiation Package  20171217
subgraph_jac_rev.hpp
Go to the documentation of this file.
3
4 /* --------------------------------------------------------------------------
6
8 the terms of the
9  Eclipse Public License Version 1.0.
10
11 A copy of this license is included in the COPYING file of this distribution.
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 Purpose$$42 We use latex F : \B{R}^n \rightarrow \B{R}^m$$ to denote the
43 function corresponding to $icode f$$. 44 Here icode n$$ is the$cref/domain/seq_property/Domain/$$size, 45 and icode m$$ is the $cref/range/seq_property/Range/$$size, or icode f$$. 46 The syntax above takes advantage of sparsity when computing the Jacobian 47$latex $48 J(x) = F^{(1)} (x) 49$ $$50 The first syntax computes the sparsity pattern and the value 51 of the Jacobian at the same time. 52 If one only wants the sparsity pattern, 53 it should be faster to use cref subgraph_sparsity$$.
54
55 $head Method$$56 This routine uses a subgraph technique. To be specific, 57 for each dependent variable, 58 it a subgraph of the operation sequence 59 to determine which independent variables affect it. 60 This avoids to overhead of performing set operations 61 that is inherent in other methods for computing sparsity patterns. 62 63 head BaseVector$$ 64 The type$icode BaseVector$$is a cref SimpleVector$$ class with
65 $cref/elements of type/SimpleVector/Elements of Specified Type/$$66 icode Base$$. 67 68$head SizeVector$$69 The type icode SizeVector$$ is a $cref SimpleVector$$class with 70 cref/elements of type/SimpleVector/Elements of Specified Type/$$ 71$code size_t$$. 72 73 head BoolVector$$
74 The type $icode BoolVector$$is a cref SimpleVector$$ class with 75$cref/elements of type/SimpleVector/Elements of Specified Type/$$76 code bool$$.
77
78 $head f$$79 This object has prototype 80 codei% 81 ADFun<%Base%> %f% 82 %$$ 83 Note that the Taylor coefficients stored in$icode f$$are affected 84 by this operation; see 85 cref/uses forward/sparse_jac/Uses Forward/$$ below.
86
87 $head x$$88 This argument has prototype 89 codei% 90 const %BaseVector%& %x% 91 %$$ 92 It is the value of$icode x$$at which we are computing the Jacobian. 93 94 head Uses Forward$$
95 After each call to $cref Forward$$, 96 the object icode f$$ contains the corresponding 97$cref/Taylor coefficients/glossary/Taylor Coefficient/$$. 98 After a call to code sparse_jac_forward$$ or $code sparse_jac_rev$$, 99 the zero order coefficients correspond to 100 codei% 101 %f%.Forward(0, %x%) 102 %$$ 103 All the other forward mode coefficients are unspecified. 104 105$head subset$$106 This argument has prototype 107 codei% 108 sparse_rcv<%SizeVector%, %BaseVector%>& %subset% 109 %$$
110 Its row size is $icode%subset%.nr() == %m%$$, 111 and its column size is icode%subset%.nc() == %n%$$. 112 It specifies which elements of the Jacobian are computed. 113 The input elements in its value vector 114$icode%subset%.val()%$$do not matter. 115 Upon return it contains the value of the corresponding elements 116 of the Jacobian. 117 118 head select_domain$$
119 The argument $icode select_domain$$has prototype 120 codei% 121 const %BoolVector%& %select_domain% 122 %$$ 123 It has size$latex n$$and specifies which independent variables 124 to include. 125 126 head select_range$$
127 The argument $icode select_range$$has prototype 128 codei% 129 const %BoolVector%& %select_range% 130 %$$ 131 It has size$latex m$$and specifies which components of the range 132 to include in the calculation. 133 A subgraph is built for each dependent variable and the selected set 134 of independent variables. 135 136 head matrix_out$$
137 This argument has prototype
138 $codei% 139 sparse_rcv<%SizeVector%, %BaseVector%>& %matrix_out% 140 %$$141 This input value of icode matrix_out$$ does not matter. 142 Upon return$icode matrix_out$$is 143 cref/sparse matrix/sparse_rcv/$$ representation of $latex F^{(1)} (x)$$. 144 The matrix has latex m$$ rows,$latex n$$columns. 145 If icode%select_domain%[%j%]%$$ is true,
146 $icode%select_range%[%i%]%$$is true, and 147 latex F_i (x)$$ depends on$latex x_j$$, 148 then the pair latex (i, j)$$ is in $icode matrix_out$$. 149 For each icode%k% = 0 , %...%, %matrix_out%.nnz()%$$, let 150$codei%
151  %i% = %matrix_out%.row()[%k%]
152  %j% = %matrix_out%.col()[%k%]
153  %v% = %matrix_out%.val()[%k%]
154 %$$155 It follows that the partial of latex F_i (x)$$ with respect to
156 $latex x_j$$is equal to latex v$$. 157 158 159$head Example$$160 children% 161 example/sparse/subgraph_jac_rev.cpp% 162 example/sparse/subgraph_hes2jac.cpp 163 %$$
164 The files $cref subgraph_jac_rev.cpp$$and cref subgraph_hes2jac.cpp$$ 165 are examples and tests using$code subgraph_jac_rev$$. 166 They returns code true$$ for success and $code false$$for failure. 167 168$end
169 -----------------------------------------------------------------------------
170 */
173
175
176 /*!
177 Subgraph sparsity patterns.
178
179 \tparam Base
180 is the base type for this recording.
181
182 \tparam SizeVector
183 is the simple vector with elements of type size_t that is used for
184 row, column index sparsity patterns.
185
186 \tparam BaseVector
187 a simple vector class with elements of type Base.
188
189 \param x
190 a vector of length n, the number of independent variables in f
192
193 \param subset
194 specifices the subset of the sparsity pattern where the Jacobian is evaluated.
195 subset.nr() == m,
196 subset.nc() == n.
197 */
198 template <typename Base>
199 template <typename SizeVector, typename BaseVector>
201  const BaseVector& x ,
203 { size_t m = Range();
204  size_t n = Domain();
205  //
207  subset.nr() == m,
208  "subgraph_jac_rev: subset.nr() not equal range dimension for f"
209  );
211  subset.nc() == n,
212  "subgraph_jac_rev: subset.nc() not equal domain dimension for f"
213  );
214  //
215  // point at which we are evaluating Jacobian
216  Forward(0, x);
217  //
218  // nnz and row, column, and row_major vectors for subset
219  size_t nnz = subset.nnz();
220  const SizeVector& row( subset.row() );
221  const SizeVector& col( subset.col() );
222  SizeVector row_major = subset.row_major();
223  //
224  // determine set of independent variabels
225  local::pod_vector<bool> select_domain(n);
226  for(size_t j = 0; j < n; j++)
227  select_domain[j] = false;
228  for(size_t k = 0; k < nnz; k++)
229  select_domain[ col[k] ] = true;
230  //
231  // initialize reverse mode computation on subgraphs
232  subgraph_reverse(select_domain);
233  //
234  // memory used to hold subgraph_reverse results
235  BaseVector dw;
236  SizeVector dw_col;
237  //
238  // initialize index in row_major
239  size_t k = 0;
240  Base zero(0);
241  while(k < nnz )
242  { size_t q = 1;
243  size_t i_dep = row[ row_major[k] ];
244  size_t i_ind = col[ row_major[k] ];
245  size_t ell = i_dep;
246  subgraph_reverse(q, ell, dw_col, dw);
247  //
248  size_t c = 0;
249  while( i_dep == ell )
250  { // check that subgraph_reverse has retured this value
251  if( c < size_t( dw_col.size() ) )
252  { if( i_ind == dw_col[c++] )
253  subset.set( row_major[k], dw[i_ind] );
254  else
255  subset.set( row_major[k], zero);
256  }
257  else
258  subset.set( row_major[k], zero);
259  ++k;
260  if( k == nnz )
261  { i_dep = m;
262  i_ind = n;
263  }
264  else
265  { i_dep = row[ row_major[k] ];
266  i_ind = col[ row_major[k] ];
267  }
268  }
269  }
270  return;
271 }
272 template <typename Base>
273 template <typename BoolVector, typename SizeVector, typename BaseVector>
275  const BoolVector& select_domain ,
276  const BoolVector& select_range ,
277  const BaseVector& x ,
279 { size_t m = Range();
280  size_t n = Domain();
281  //
282  // point at which we are evaluating Jacobian
283  Forward(0, x);
284  //
285  // nnz and row, column, and row_major vectors for subset
288  local::pod_vector<Base> val_out;
289  //
290  // initialize reverse mode computation on subgraphs
291  subgraph_reverse(select_domain);
292  //
293  // memory used to hold subgraph_reverse results
294  BaseVector dw;
295  SizeVector col;
296  //
297  // loop through selected independent variables
298  for(size_t i = 0; i < m; ++i) if( select_range[i] )
299  { // compute Jacobian and sparsity for this dependent variable
300  size_t q = 1;
301  subgraph_reverse(q, i, col, dw);
302  CPPAD_ASSERT_UNKNOWN( size_t( dw.size() ) == n );
303  //
304  // offset for this dependent variable
305  size_t index = row_out.size();
306  CPPAD_ASSERT_UNKNOWN( col_out.size() == index );
307  CPPAD_ASSERT_UNKNOWN( val_out.size() == index );
308  //
309  // extend vectors to hold results for this dependent variable
310  size_t col_size = size_t( col.size() );
311  row_out.extend( col_size );
312  col_out.extend( col_size );
313  val_out.extend( col_size );
314  //
315  // store results for this dependent variable
316  for(size_t c = 0; c < col_size; ++c)
317  { row_out[index + c] = i;
318  col_out[index + c] = col[c];
319  val_out[index + c] = dw[ col[c] ];
320  }
321  }
322  //
323  // create sparsity pattern corresponding to row_out, col_out
324  size_t nr = m;
325  size_t nc = n;
326  size_t nnz = row_out.size();
327  sparse_rc<SizeVector> pattern(nr, nc, nnz);
328  for(size_t k = 0; k < nnz; ++k)
329  pattern.set(k, row_out[k], col_out[k]);
330  //
331  // create sparse matrix
332  sparse_rcv<SizeVector, BaseVector> matrix(pattern);
333  for(size_t k = 0; k < nnz; ++k)
334  matrix.set(k, val_out[k]);
335  //
336  // return matrix
337  matrix_out = matrix;
338  //
339  return;
340 }
342 # endif
Check that exp is true, if not print msg and terminate execution.
size_t extend(size_t n)
Increase the number of elements the end of this vector (existing elements are always preserved)...
Definition: pod_vector.hpp:115
void set(size_t k, const value_type &v)
Definition: sparse_rcv.hpp:240
subgraph information attached to a operation sequence
size_t nnz(void) const
number of possibly non-zero elements in matrix
Definition: sparse_rcv.hpp:254
const SizeVector & col(void) const
column indices
Definition: sparse_rcv.hpp:260
void subgraph_jac_rev(const BaseVector &x, sparse_rcv< SizeVector, BaseVector > &subset)
Subgraph sparsity patterns.
size_t nc(void) const
number of columns in matrix
Definition: sparse_rcv.hpp:251
Sparse matrices with elements of type Scalar.
Definition: sparse_rcv.hpp:212
size_t nr(void) const
number of rows in matrix
Definition: sparse_rcv.hpp:248
sparsity pattern for a matrix with indices of type size_t
Definition: sparse_rc.hpp:212
size_t size(void) const
current number of elements in this vector.
Definition: pod_vector.hpp:79
SizeVector row_major(void) const
row-major order
Definition: sparse_rcv.hpp:266