1 | /* $Id: rev_sparse_jac.hpp 1532 2009-09-28 11:55:30Z bradbell $ */ |
---|
2 | # ifndef CPPAD_REV_SPARSE_JAC_INCLUDED |
---|
3 | # define CPPAD_REV_SPARSE_JAC_INCLUDED |
---|
4 | |
---|
5 | /* -------------------------------------------------------------------------- |
---|
6 | CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 Bradley M. Bell |
---|
7 | |
---|
8 | CppAD is distributed under multiple licenses. This distribution is under |
---|
9 | the terms of the |
---|
10 | Common Public License Version 1.0. |
---|
11 | |
---|
12 | A copy of this license is included in the COPYING file of this distribution. |
---|
13 | Please visit http://www.coin-or.org/CppAD/ for information on other licenses. |
---|
14 | -------------------------------------------------------------------------- */ |
---|
15 | |
---|
16 | /* |
---|
17 | $begin RevSparseJac$$ |
---|
18 | $spell |
---|
19 | VecAD |
---|
20 | var |
---|
21 | Jacobian |
---|
22 | Jac |
---|
23 | const |
---|
24 | Bool |
---|
25 | Dep |
---|
26 | proportional |
---|
27 | $$ |
---|
28 | |
---|
29 | $section Jacobian Sparsity Pattern: Reverse Mode$$ |
---|
30 | |
---|
31 | $index RevSparseJac$$ |
---|
32 | $index reverse, sparse Jacobian$$ |
---|
33 | $index sparse, reverse Jacobian$$ |
---|
34 | $index pattern, reverse Jacobian$$ |
---|
35 | |
---|
36 | $head Syntax$$ |
---|
37 | $syntax%%r% = %F%.RevSparseJac(%p%, %s%)%$$ |
---|
38 | $pre |
---|
39 | $$ |
---|
40 | $syntax%%r% = %F%.RevSparseJac(%p%, %s%, %packed%)%$$ |
---|
41 | |
---|
42 | |
---|
43 | $head Purpose$$ |
---|
44 | We use $latex F : B^n \rightarrow B^m$$ to denote the |
---|
45 | $xref/glossary/AD Function/AD function/$$ corresponding to $italic f$$. |
---|
46 | For a fixed $latex p \times m$$ matrix $latex S$$, |
---|
47 | the Jacobian of $latex S * F( x )$$ |
---|
48 | with respect to $latex x$$ is |
---|
49 | $latex \[ |
---|
50 | J(x) = S * F^{(1)} ( x ) |
---|
51 | \] $$ |
---|
52 | Given a |
---|
53 | $xref/glossary/Sparsity Pattern/sparsity pattern/$$ |
---|
54 | for $latex S$$, |
---|
55 | $code RevSparseJac$$ returns a sparsity pattern for the $latex J(x)$$. |
---|
56 | |
---|
57 | $head f$$ |
---|
58 | The object $italic f$$ has prototype |
---|
59 | $syntax% |
---|
60 | ADFun<%Base%> %f% |
---|
61 | %$$ |
---|
62 | |
---|
63 | $head x$$ |
---|
64 | the sparsity pattern is valid for all values of the independent |
---|
65 | variables in $latex x \in B^n$$ |
---|
66 | (even if you use $cref/CondExp/$$ or $cref/VecAD/$$). |
---|
67 | |
---|
68 | $head p$$ |
---|
69 | The argument $italic p$$ has prototype |
---|
70 | $syntax% |
---|
71 | size_t %p% |
---|
72 | %$$ |
---|
73 | It specifies the number of rows in the Jacobian $latex J(x)$$. |
---|
74 | Note that the memory required for the calculation is proportional |
---|
75 | to $latex p$$ times the total number of variables |
---|
76 | in the AD operation sequence corresponding to $italic f$$ |
---|
77 | ($xref/SeqProperty/size_var/f.size_var/$$). |
---|
78 | Smaller values for $italic p$$ can be used to |
---|
79 | break the sparsity calculation |
---|
80 | into groups that do not require to much memory. |
---|
81 | |
---|
82 | $head s$$ |
---|
83 | The argument $italic s$$ has prototype |
---|
84 | $syntax% |
---|
85 | const %VectorBool% &%s% |
---|
86 | %$$ |
---|
87 | (see $xref/RevSparseJac/VectorBool/VectorBool/$$ below) |
---|
88 | and its size is $latex p * m$$. |
---|
89 | It specifies a |
---|
90 | $xref/glossary/Sparsity Pattern/sparsity pattern/$$ |
---|
91 | for the matrix $italic S$$ as follows: |
---|
92 | for $latex i = 0 , \ldots , p-1$$ and $latex j = 0 , \ldots , m-1$$. |
---|
93 | $latex \[ |
---|
94 | S_{i,j} \neq 0 ; \Rightarrow \; s [ i * m + j ] = {\rm true} |
---|
95 | \] $$ |
---|
96 | |
---|
97 | $head packed$$ |
---|
98 | If $italic packed$$ is true, |
---|
99 | during the sparsity calculation sets of indices are represented |
---|
100 | as vectors of bits that packed into words and operations are done |
---|
101 | on multiple bits at a time (the number of bits in a word is unspecified). |
---|
102 | Otherwise, sets of indices are represents using a sparse structure |
---|
103 | that only includes the non-zero indices and operations are done |
---|
104 | one index at a time. |
---|
105 | $pre |
---|
106 | |
---|
107 | $$ |
---|
108 | The default value for $italic packed$$ is true; i.e., |
---|
109 | the value used if it is not present. |
---|
110 | |
---|
111 | $head r$$ |
---|
112 | The return value $italic r$$ has prototype |
---|
113 | $syntax% |
---|
114 | %VectorBool% %r% |
---|
115 | %$$ |
---|
116 | (see $xref/RevSparseJac/VectorBool/VectorBool/$$ below) |
---|
117 | and its size is $latex p * n$$. |
---|
118 | It specifies a |
---|
119 | $xref/glossary/Sparsity Pattern/sparsity pattern/$$ |
---|
120 | for the matrix $latex J(x)$$ as follows: |
---|
121 | for $latex x \in B^n$$, |
---|
122 | for $latex i = 0 , \ldots , p-1$$, |
---|
123 | and $latex j = 0 , \ldots , n-1$$ |
---|
124 | $latex \[ |
---|
125 | J(x)_{i,j} \neq 0 ; \Rightarrow \; r [ i * n + j ] = {\rm true} |
---|
126 | \] $$ |
---|
127 | |
---|
128 | $head VectorBool$$ |
---|
129 | The type $italic VectorBool$$ must be a $xref/SimpleVector/$$ class with |
---|
130 | $xref/SimpleVector/Elements of Specified Type/elements of type bool/$$. |
---|
131 | The routine $xref/CheckSimpleVector/$$ will generate an error message |
---|
132 | if this is not the case. |
---|
133 | In order to save memory, |
---|
134 | you may want to use a class that packs multiple elements into one |
---|
135 | storage location; for example, |
---|
136 | $xref/CppAD_vector/vectorBool/vectorBool/$$. |
---|
137 | |
---|
138 | |
---|
139 | $head Entire Sparsity Pattern$$ |
---|
140 | Suppose that $latex p = m$$ and |
---|
141 | $latex S$$ is the $latex m \times m$$ identity matrix, |
---|
142 | If follows that |
---|
143 | $latex \[ |
---|
144 | s [ i * q + j ] = \left\{ \begin{array}{ll} |
---|
145 | {\rm true} & {\rm if} \; i = j \\ |
---|
146 | {\rm flase} & {\rm otherwise} |
---|
147 | \end{array} \right. |
---|
148 | \] $$ |
---|
149 | is an efficient sparsity pattern for the Jacobian of $latex S$$; |
---|
150 | i.e., the choice for $italic s$$ has as few true values as possible. |
---|
151 | In this case, |
---|
152 | the corresponding value for $italic r$$ is a |
---|
153 | sparsity pattern for the Jacobian $latex J(x) = F^{(1)} ( x )$$. |
---|
154 | |
---|
155 | |
---|
156 | $head Example$$ |
---|
157 | $children% |
---|
158 | example/rev_sparse_jac.cpp |
---|
159 | %$$ |
---|
160 | The file |
---|
161 | $xref/RevSparseJac.cpp/$$ |
---|
162 | contains an example and test of this operation. |
---|
163 | It returns true if it succeeds and false otherwise. |
---|
164 | |
---|
165 | $end |
---|
166 | ----------------------------------------------------------------------------- |
---|
167 | */ |
---|
168 | |
---|
169 | // BEGIN CppAD namespace |
---|
170 | namespace CppAD { |
---|
171 | |
---|
172 | template <class Base, class VectorBool, class VectorSet> |
---|
173 | void ForSparseJac( |
---|
174 | size_t total_num_var , |
---|
175 | size_t p , |
---|
176 | const VectorBool& s , |
---|
177 | CppAD::vector<size_t>& dep_taddr , |
---|
178 | CppAD::vector<size_t>& ind_taddr , |
---|
179 | CppAD::player<Base>& play , |
---|
180 | VectorBool& r ) |
---|
181 | { |
---|
182 | // temporary indices |
---|
183 | size_t i, j; |
---|
184 | |
---|
185 | // check VectorBool is Simple Vector class with bool elements |
---|
186 | CheckSimpleVector<bool, VectorBool>(); |
---|
187 | |
---|
188 | // range and domain dimensions for F |
---|
189 | size_t m = dep_taddr.size(); |
---|
190 | size_t n = ind_taddr.size(); |
---|
191 | |
---|
192 | CPPAD_ASSERT_KNOWN( |
---|
193 | p > 0, |
---|
194 | "RevSparseJac: p (first argument) is not greater than zero" |
---|
195 | ); |
---|
196 | |
---|
197 | CPPAD_ASSERT_KNOWN( |
---|
198 | s.size() == p * m, |
---|
199 | "RevSparseJac: s (second argument) length is not equal to\n" |
---|
200 | "p (first argument) times range dimension for ADFun object." |
---|
201 | ); |
---|
202 | |
---|
203 | // vector of sets that will hold the results |
---|
204 | VectorSet var_sparsity; |
---|
205 | var_sparsity.resize(total_num_var, p); |
---|
206 | |
---|
207 | // The sparsity pattern corresponding to the dependent variables |
---|
208 | for(i = 0; i < m; i++) |
---|
209 | { CPPAD_ASSERT_UNKNOWN( dep_taddr[i] < total_num_var ); |
---|
210 | |
---|
211 | for(j = 0; j < p; j++) if( s[ i * m + j ] ) |
---|
212 | var_sparsity.add_element( dep_taddr[i], j ); |
---|
213 | } |
---|
214 | |
---|
215 | // evaluate the sparsity patterns |
---|
216 | RevJacSweep( |
---|
217 | n, |
---|
218 | total_num_var, |
---|
219 | &play, |
---|
220 | var_sparsity |
---|
221 | ); |
---|
222 | |
---|
223 | // return values corresponding to dependent variables |
---|
224 | CPPAD_ASSERT_UNKNOWN( r.size() == p * n ); |
---|
225 | for(j = 0; j < n; j++) |
---|
226 | { CPPAD_ASSERT_UNKNOWN( ind_taddr[j] == (j+1) ); |
---|
227 | |
---|
228 | // ind_taddr[j] is operator taddr for j-th independent variable |
---|
229 | CPPAD_ASSERT_UNKNOWN( play.GetOp( ind_taddr[j] ) == InvOp ); |
---|
230 | |
---|
231 | // set bits |
---|
232 | for(i = 0; i < p; i++) |
---|
233 | r[ i * n + j ] = false; |
---|
234 | |
---|
235 | CPPAD_ASSERT_UNKNOWN( var_sparsity.limit() == p ); |
---|
236 | i = var_sparsity.next_element(j+1); |
---|
237 | while( i < p ) |
---|
238 | { r[ i * n + j ] = true; |
---|
239 | i = var_sparsity.next_element(j+1); |
---|
240 | } |
---|
241 | } |
---|
242 | } |
---|
243 | |
---|
244 | template <class Base> |
---|
245 | template <class VectorBool> |
---|
246 | VectorBool ADFun<Base>::RevSparseJac( |
---|
247 | size_t p , |
---|
248 | const VectorBool& s , |
---|
249 | bool packed ) |
---|
250 | { size_t n = ind_taddr_.size(); |
---|
251 | VectorBool r( p * n ); |
---|
252 | |
---|
253 | if( packed ) |
---|
254 | { CppAD::ForSparseJac<Base, VectorBool, vector_pack>( |
---|
255 | total_num_var_ , |
---|
256 | p , |
---|
257 | s , |
---|
258 | dep_taddr_ , |
---|
259 | ind_taddr_ , |
---|
260 | play_ , |
---|
261 | r |
---|
262 | ); |
---|
263 | } |
---|
264 | else |
---|
265 | { CppAD::ForSparseJac<Base, VectorBool, vector_set>( |
---|
266 | total_num_var_ , |
---|
267 | p , |
---|
268 | s , |
---|
269 | dep_taddr_ , |
---|
270 | ind_taddr_ , |
---|
271 | play_ , |
---|
272 | r |
---|
273 | ); |
---|
274 | } |
---|
275 | return r; |
---|
276 | } |
---|
277 | |
---|
278 | } // END CppAD namespace |
---|
279 | |
---|
280 | |
---|
281 | # endif |
---|