source: trunk/cppad/core/sparse_hessian.hpp @ 3941

Last change on this file since 3941 was 3941, checked in by bradbell, 2 years ago

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

commit c8c4cc081accff3628e7e66370ec01e4c99afe8d
Author: Brad Bell <bradbell@…>
Date: Thu Jun 1 23:16:39 2017 -0600

Changes automatically generated by the autotools.

commit f4392bc3eee8f6d0ccd45a0bb3be51181e211680
Author: Brad Bell <bradbell@…>
Date: Thu Jun 1 23:11:56 2017 -0600

  1. Add colpack_jac.cpp example (rename colpack_jac.cpp->colpack_jacobian.cpp).
  2. Add colpack_hescpp example (rename colpack_hes.cpp->colpack_hessian.cpp).


test_one.sh.in: adapt to using test_boolofvoid for testing.
sparse_hes.hpp: fix bug in cppad.symmetric case.

commit 086b8a8709b0c9cb01ce2cf8bc7910e903105ff7
Author: Brad Bell <bradbell@…>
Date: Thu Jun 1 08:54:59 2017 -0600

  1. Fix bug in use of colpack (see kludge in comments).
  2. Fix colpack.symmetric (not general) and add colpack.general.
  3. Deprecate colpack.star.
  4. More autotools from install to deprecated.
  5. Advance to cppad-20170601.

commit 23f26c060648f5c6fc62a1598c659aeccc5ca46f
Author: Brad Bell <bradbell@…>
Date: Tue May 30 08:14:04 2017 -0700

Advance to cppad-20170530.

commit 97f8c08509865d1bfb7ec2e5cd557ddc979f8412
Author: Brad Bell <bradbell@…>
Date: Tue May 30 07:38:47 2017 -0700

debug_rel branch:
There is a problem with speed sparse_hessian debug that goes back to master.
Supresss debug in cppad speed tests until it is fixed.

commit 39ea0d7d9c041784ccd26ce80d19a7ab02752818
Author: Brad Bell <bradbell@…>
Date: Mon May 29 22:34:22 2017 -0700

debug_rel branch:
run_cmake.sh: fix debug_none case.
CMakeLists.txt: use cppad_debug_which to determine debug or release.
CMakeLists.txt: let set_compile_flags determkine build type.

commit 191553e54dca407207789cf0d7c6c27fe6188775
Author: Brad Bell <bradbell@…>
Date: Mon May 29 19:53:08 2017 -0700

debug_rel branch:
Use set_compile_flags in introduction.

commit fba276a84e58d9a0d0944168d5706b7446beb32c
Author: Brad Bell <bradbell@…>
Date: Mon May 29 19:46:30 2017 -0700

debug_rel branch:
Use set_compile_flags in eample/multi_thread subdirectories.

commit 66c8cdb266fa3af29b211b8c870a3aed7a13b021
Author: Brad Bell <bradbell@…>
Date: Mon May 29 18:56:48 2017 -0700

debug_rel branch:
Use set_compile_flags in speed directory.

commit c431b15ee7714d3106234bc527ba2f9a836750e1
Author: Brad Bell <bradbell@…>
Date: Mon May 29 18:36:51 2017 -0700

debug_rel branch:
Convert cppad_ipopt to use set_compile_flags and cppad_debug_which.


CMakeLists.txt: alwasy compile for release to reduce testing time.

commit 2c95b0019f1b665fb14b9f00b049e8b5fb11f89d
Author: Brad Bell <bradbell@…>
Date: Mon May 29 16:55:07 2017 -0700

debug_rel branch:
Add cppad_debug_which to the cmake command line.

commit fd8d1498cf6dc092deca41f764cbb2a001a4cc88
Author: Brad Bell <bradbell@…>
Date: Mon May 29 08:14:23 2017 -0700

debug_rel branch:
Change random_debug_release -> set_compile_flags.

commit 159f5a5aa09012213a52f4ed1c9f0607129a5fe7
Author: Brad Bell <bradbell@…>
Date: Mon May 29 06:50:43 2017 -0700

debug_rel branch:
Update the autotools automatically generated build files.


batch_edit.sh: Start comments about a plan for editing all the source files.
get_sacado.sh: advance to trilions-12.10.11.
makefile.am: advance to trilinos-12.10.1

commit 302153317cd296ec6f927c3202cf96bf38594bbb
Author: Brad Bell <bradbell@…>
Date: Mon May 29 05:20:00 2017 -0700

debug_rel branch:
Add error message if sacado configuration file does not exist.

commit 3f01a631ae808c3a1359e53e1cd55e9a0ea88711
Author: Brad Bell <bradbell@…>
Date: Mon May 29 04:24:00 2017 -0700

debug_rel branch:
CMakeLists.txt: automate naming of libraries Sacado needs.
checkpoint.cpp: fix warnings.

commit dd240928c0c8b6972a8197c985ccc01f08b8886b
Author: Brad Bell <bradbell@…>
Date: Sun May 28 08:25:20 2017 -0700

debug_rel branch:
sparse_sub_hes.cpp: add missing cases found by clang compiler.

commit 30a0c35f1ac50ec425be9a2b7b026284026eccd7
Author: Brad Bell <bradbell@…>
Date: Sun May 28 07:57:36 2017 -0700

debug_rel branch:
eigen_cholesky.hpp: fix compiler warning.
harmonic_time.cpp: remove include that is not used.
forward_active.cpp: fix compiler warning.

commit 4876d14e49dc235865b1574fb38a55cf5ea7a422
Author: Brad Bell <bradbell@…>
Date: Sun May 28 06:19:48 2017 -0700

debug_rel branch:
random_debug_release.cmake: fix comment, remove message replaced by random_choice_0123 in output.
multiple_solution.cpp: fix warnings with clang compiler.
eigen_cholesky.hpp: fix warnings with clang compiler.
compare_change.cpp: fix CPPAD_DEBUG_AND_RELEASE case.

commit 2c51a18f35188d04d2f94069382439580e23f4ac
Author: Brad Bell <bradbell@…>
Date: Sat May 27 21:04:37 2017 -0700

debug_rel branch:
Advance version to cppad-20170527.

commit 4500887b362537774b05e954ad2a95b65a7b8ba0
Author: Brad Bell <bradbell@…>
Date: Sat May 27 09:04:56 2017 -0700

debug_rel branch:
Ramdomly select debug or release flags in example directory.


CMakeLists.txt: always debug for multi_threed examples.

commit 140b5269a0b1a30643894e5a7a8c9a5eb1310301
Author: Brad Bell <bradbell@…>
Date: Sat May 27 08:10:51 2017 -0700

debug_rel branch:
Changing how we set all debug and release flags.

commit e6fb2639db1288fb75de4030b5906df1e41756f9
Author: Brad Bell <bradbell@…>
Date: Sat May 27 07:30:24 2017 -0700

debug_rel branch:
Replace use of cppad_extra_debug by CPPAD_DEBUG_AND_RELEASE.

commit fbbfd0f6e94862174a8a7a17308489ffddb28084
Author: Brad Bell <bradbell@…>
Date: Sat May 27 05:55:58 2017 -0700

debug_rel branch:
Improve random selection of which files are build for release or debug.


forward.cpp: use new -DCPPAD_DEBUG_AND_RELEASE flag.

commit 284be366fb5e2f685a0c71ea6a0e3f74584bf187
Author: Brad Bell <bradbell@…>
Date: Thu May 25 07:39:32 2017 -0700

debug_rel branch:
Add test that failed before change to player.


player.hpp: Fix so it has the same size in debug and release more.
checkpoint.cpp: fix warning when compiling for release.
run_cmake.sh: prepare to use random number to switch debug and release set.
CMakeLists.txt: switch to only test debug (for now).

commit f32375b77e3825628fee6cb160f691a32c48b796
Author: Brad Bell <bradbell@…>
Date: Wed May 24 12:04:27 2017 -0700

debug_rel branch:
forward.cpp: fix a warning during release build.

commit 5fcc7eb78ae8de9f1dbc6c4f0c76fe38e8aeba95
Author: Brad Bell <bradbell@…>
Date: Wed May 24 10:11:12 2017 -0700

debug_rel branch:
CMakeLists.txt: make easy to mix debug and release builds.
eigen_mat_inv.hpp: fix release version warning.

commit 696266f3d62079f5e3bfb1a0f60a7e4f8134e068
Author: Brad Bell <bradbell@…>
Date: Wed May 24 05:43:29 2017 -0700

push_git2svn.py: user ./build in place of ./build/work.
testvector.hpp: improve comments about replacing CPPAD_TESTVECTOR.

File size: 25.0 KB
Line 
1# ifndef CPPAD_CORE_SPARSE_HESSIAN_HPP
2# define CPPAD_CORE_SPARSE_HESSIAN_HPP
3
4/* --------------------------------------------------------------------------
5CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
6
7CppAD is distributed under multiple licenses. This distribution is under
8the terms of the
9                    Eclipse Public License Version 1.0.
10
11A copy of this license is included in the COPYING file of this distribution.
12Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
13-------------------------------------------------------------------------- */
14
15/*
16$begin sparse_hessian$$
17$spell
18        jacobian
19        recomputed
20        CppAD
21        valarray
22        std
23        Bool
24        hes
25        const
26        Taylor
27        cppad
28        cmake
29        colpack
30$$
31
32$section Sparse Hessian$$
33
34$head Syntax$$
35$icode%hes% = %f%.SparseHessian(%x%, %w%)
36%hes% = %f%.SparseHessian(%x%, %w%, %p%)
37%n_sweep% = %f%.SparseHessian(%x%, %w%, %p%, %row%, %col%, %hes%, %work%)
38%$$
39
40$head Purpose$$
41We use $latex n$$ for the $cref/domain/seq_property/Domain/$$ size,
42and $latex m$$ for the $cref/range/seq_property/Range/$$ size of $icode f$$.
43We use $latex F : \B{R}^n \rightarrow \B{R}^m$$ do denote the
44$cref/AD function/glossary/AD Function/$$
45corresponding to $icode f$$.
46The syntax above sets $icode hes$$ to the Hessian
47$latex \[
48        H(x) = \dpow{2}{x} \sum_{i=1}^m w_i F_i (x)
49\] $$
50This routine takes advantage of the sparsity of the Hessian
51in order to reduce the amount of computation necessary.
52If $icode row$$ and $icode col$$ are present, it also takes
53advantage of the reduced set of elements of the Hessian that
54need to be computed.
55One can use speed tests (e.g. $cref speed_test$$)
56to verify that results are computed faster
57than when using the routine $cref Hessian$$.
58
59$head f$$
60The object $icode f$$ has prototype
61$codei%
62        ADFun<%Base%> %f%
63%$$
64Note that the $cref ADFun$$ object $icode f$$ is not $code const$$
65(see $cref/Uses Forward/sparse_hessian/Uses Forward/$$ below).
66
67$head x$$
68The argument $icode x$$ has prototype
69$codei%
70        const %VectorBase%& %x%
71%$$
72(see $cref/VectorBase/sparse_hessian/VectorBase/$$ below)
73and its size
74must be equal to $icode n$$, the dimension of the
75$cref/domain/seq_property/Domain/$$ space for $icode f$$.
76It specifies
77that point at which to evaluate the Hessian.
78
79$head w$$
80The argument $icode w$$ has prototype
81$codei%
82        const %VectorBase%& %w%
83%$$
84and size $latex m$$.
85It specifies the value of $latex w_i$$ in the expression
86for $icode hes$$.
87The more components of $latex w$$ that are identically zero,
88the more sparse the resulting Hessian may be (and hence the more efficient
89the calculation of $icode hes$$ may be).
90
91$head p$$
92The argument $icode p$$ is optional and has prototype
93$codei%
94        const %VectorSet%& %p%
95%$$
96(see $cref/VectorSet/sparse_hessian/VectorSet/$$ below)
97If it has elements of type $code bool$$,
98its size is $latex n * n$$.
99If it has elements of type $code std::set<size_t>$$,
100its size is $latex n$$ and all its set elements are between
101zero and $latex n - 1$$.
102It specifies a
103$cref/sparsity pattern/glossary/Sparsity Pattern/$$
104for the Hessian $latex H(x)$$.
105
106$subhead Purpose$$
107If this sparsity pattern does not change between calls to
108$codei SparseHessian$$, it should be faster to calculate $icode p$$ once and
109pass this argument to $codei SparseHessian$$.
110If you specify $icode p$$, CppAD will use the same
111type of sparsity representation
112(vectors of $code bool$$ or vectors of $code std::set<size_t>$$)
113for its internal calculations.
114Otherwise, the representation
115for the internal calculations is unspecified.
116
117$subhead work$$
118If you specify $icode work$$ in the calling sequence,
119it is not necessary to keep the sparsity pattern; see the heading
120$cref/p/sparse_hessian/work/p/$$ under the $icode work$$ description.
121
122$subhead Column Subset$$
123If the arguments $icode row$$ and $icode col$$ are present,
124and $cref/color_method/sparse_hessian/work/color_method/$$ is
125$code cppad.general$$ or $code cppad.symmetric$$,
126it is not necessary to compute the entire sparsity pattern.
127Only the following subset of column values will matter:
128$codei%
129        { %col%[%k%] : %k% = 0 , %...% , %K%-1 }
130%$$.
131
132
133$head row, col$$
134The arguments $icode row$$ and $icode col$$ are optional and have prototype
135$codei%
136        const %VectorSize%& %row%
137        const %VectorSize%& %col%
138%$$
139(see $cref/VectorSize/sparse_hessian/VectorSize/$$ below).
140They specify which rows and columns of $latex H (x)$$ are
141returned and in what order.
142We use $latex K$$ to denote the value $icode%hes%.size()%$$
143which must also equal the size of $icode row$$ and $icode col$$.
144Furthermore,
145for $latex k = 0 , \ldots , K-1$$, it must hold that
146$latex row[k] < n$$ and $latex col[k] < n$$.
147In addition,
148all of the $latex (row[k], col[k])$$ pairs must correspond to a true value
149in the sparsity pattern $icode p$$.
150
151$head hes$$
152The result $icode hes$$ has prototype
153$codei%
154        %VectorBase% %hes%
155%$$
156In the case where $icode row$$ and $icode col$$ are not present,
157the size of $icode hes$$ is $latex n * n$$ and
158its size is $latex n * n$$.
159In this case, for $latex i = 0 , \ldots , n - 1 $$
160and $latex ell = 0 , \ldots , n - 1$$
161$latex \[
162        hes [ j * n + \ell ] = \DD{ w^{\rm T} F }{ x_j }{ x_\ell } ( x )
163\] $$
164$pre
165
166$$
167In the case where the arguments $icode row$$ and $icode col$$ are present,
168we use $latex K$$ to denote the size of $icode hes$$.
169The input value of its elements does not matter.
170Upon return, for $latex k = 0 , \ldots , K - 1$$,
171$latex \[
172        hes [ k ] = \DD{ w^{\rm T} F }{ x_j }{ x_\ell } (x)
173        \; , \;
174        \; {\rm where} \;
175        j = row[k]
176        \; {\rm and } \;
177        \ell = col[k]
178\] $$
179
180$head work$$
181If this argument is present, it has prototype
182$codei%
183        sparse_hessian_work& %work%
184%$$
185This object can only be used with the routines $code SparseHessian$$.
186During its the first use, information is stored in $icode work$$.
187This is used to reduce the work done by future calls to $code SparseHessian$$
188with the same $icode f$$, $icode p$$, $icode row$$, and $icode col$$.
189If a future call is made where any of these values have changed,
190you must first call $icode%work%.clear()%$$
191to inform CppAD that this information needs to be recomputed.
192
193$subhead color_method$$
194The coloring algorithm determines which rows and columns
195can be computed during the same sweep.
196This field has prototype
197$codei%
198        std::string %work%.color_method
199%$$
200This value only matters on the first call to $code sparse_hessian$$ that
201follows the $icode work$$ constructor or a call to
202$icode%work%.clear()%$$.
203$codei%
204
205"cppad.symmetric"
206%$$
207This is the default coloring method (after a constructor or $code clear()$$).
208It takes advantage of the fact that the Hessian matrix
209is symmetric to find a coloring that requires fewer
210$cref/sweeps/sparse_hessian/n_sweep/$$.
211$codei%
212
213"cppad.general"
214%$$
215This is the same as the $code "cppad"$$ method for the
216$cref/sparse_jacobian/sparse_jacobian/work/color_method/$$ calculation.
217$codei%
218
219"colpack.symmetric"
220%$$
221This method requires that
222$cref colpack_prefix$$ was specified on the
223$cref/cmake command/cmake/CMake Command/$$ line.
224It also takes advantage of the fact that the Hessian matrix is symmetric.
225$codei%
226
227"colpack.general"
228%$$
229This is the same as the $code "colpack"$$ method for the
230$cref/sparse_jacobian/sparse_jacobian/work/color_method/$$ calculation.
231
232$subhead colpack.star Deprecated 2017-06-01$$
233The $code colpack.star$$ method is deprecated.
234It is the same as the $code colpack.symmetric$$
235which should be used instead.
236
237$subhead p$$
238If $icode work$$ is present, and it is not the first call after
239its construction or a clear,
240the sparsity pattern $icode p$$ is not used.
241This enables one to free the sparsity pattern
242and still compute corresponding sparse Hessians.
243
244$head n_sweep$$
245The return value $icode n_sweep$$ has prototype
246$codei%
247        size_t %n_sweep%
248%$$
249It is the number of first order forward sweeps
250used to compute the requested Hessian values.
251Each first forward sweep is followed by a second order reverse sweep
252so it is also the number of reverse sweeps.
253This is proportional to the total work that $code SparseHessian$$ does,
254not counting the zero order forward sweep,
255or the work to combine multiple columns into a single
256forward-reverse sweep pair.
257
258$head VectorBase$$
259The type $icode VectorBase$$ must be a $cref SimpleVector$$ class with
260$cref/elements of type/SimpleVector/Elements of Specified Type/$$
261$icode Base$$.
262The routine $cref CheckSimpleVector$$ will generate an error message
263if this is not the case.
264
265$head VectorSet$$
266The type $icode VectorSet$$ must be a $cref SimpleVector$$ class with
267$cref/elements of type/SimpleVector/Elements of Specified Type/$$
268$code bool$$ or $code std::set<size_t>$$;
269see $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for a discussion
270of the difference.
271The routine $cref CheckSimpleVector$$ will generate an error message
272if this is not the case.
273
274$subhead Restrictions$$
275If $icode VectorSet$$ has elements of $code std::set<size_t>$$,
276then $icode%p%[%i%]%$$ must return a reference (not a copy) to the
277corresponding set.
278According to section 26.3.2.3 of the 1998 C++ standard,
279$code std::valarray< std::set<size_t> >$$ does not satisfy
280this condition.
281
282$head VectorSize$$
283The type $icode VectorSize$$ must be a $cref SimpleVector$$ class with
284$cref/elements of type/SimpleVector/Elements of Specified Type/$$
285$code size_t$$.
286The routine $cref CheckSimpleVector$$ will generate an error message
287if this is not the case.
288
289$head Uses Forward$$
290After each call to $cref Forward$$,
291the object $icode f$$ contains the corresponding
292$cref/Taylor coefficients/glossary/Taylor Coefficient/$$.
293After a call to any of the sparse Hessian routines,
294the zero order Taylor coefficients correspond to
295$icode%f%.Forward(0, %x%)%$$
296and the other coefficients are unspecified.
297
298$children%
299        example/sparse/sparse_hessian.cpp%
300        example/sparse/sub_sparse_hes.cpp%
301        example/sparse/sparse_sub_hes.cpp
302%$$
303
304$head Example$$
305The routine
306$cref sparse_hessian.cpp$$
307is examples and tests of $code sparse_hessian$$.
308It return $code true$$, if it succeeds and $code false$$ otherwise.
309
310$head Subset Hessian$$
311The routine
312$cref sub_sparse_hes.cpp$$
313is an example and test that compute a sparse Hessian
314for a subset of the variables.
315It returns $code true$$, for success, and $code false$$ otherwise.
316
317$end
318-----------------------------------------------------------------------------
319*/
320# include <cppad/local/std_set.hpp>
321# include <cppad/local/color_general.hpp>
322# include <cppad/local/color_symmetric.hpp>
323
324namespace CppAD { // BEGIN_CPPAD_NAMESPACE
325/*!
326\file sparse_hessian.hpp
327Sparse Hessian driver routine and helper functions.
328*/
329// ===========================================================================
330/*!
331class used by SparseHessian to hold information
332so it does not need to be recomputed.
333*/
334class sparse_hessian_work {
335        public:
336                /// Coloring method: "cppad", or "colpack"
337                /// (this field is set by user)
338                std::string color_method;
339                /// row and column indicies for return values
340                /// (some may be reflected by star coloring algorithm)
341                CppAD::vector<size_t> row;
342                CppAD::vector<size_t> col;
343                /// indices that sort the user row and col arrays by color
344                CppAD::vector<size_t> order;
345                /// results of the coloring algorithm
346                CppAD::vector<size_t> color;
347
348                /// constructor
349                sparse_hessian_work(void) : color_method("cppad.symmetric")
350                { }
351                /// inform CppAD that this information needs to be recomputed
352                void clear(void)
353                {       color_method = "cppad.symmetric";
354                        row.clear();
355                        col.clear();
356                        order.clear();
357                        color.clear();
358                }
359};
360// ===========================================================================
361/*!
362Private helper function that does computation for all Sparse Hessian cases.
363
364\tparam Base
365is the base type for the recording that is stored in this ADFun<Base object.
366
367\tparam VectorBase
368is a simple vector class with elements of type \a Base.
369
370\tparam VectorSet
371is a simple vector class with elements of type
372\c bool or \c std::set<size_t>.
373
374\tparam VectorSize
375is sparse_pack or sparse_list.
376
377\param x [in]
378is a vector specifing the point at which to compute the Hessian.
379
380\param w [in]
381is the weighting vector that defines a scalar valued function by
382a weighted sum of the components of the vector valued function
383$latex F(x)$$.
384
385\param sparsity [in]
386is the sparsity pattern for the Hessian that we are calculating.
387
388\param user_row [in]
389is the vector of row indices for the returned Hessian values.
390
391\param user_col [in]
392is the vector of columns indices for the returned Hessian values.
393It must have the same size as user_row.
394
395\param hes [out]
396is the vector of Hessian values.
397It must have the same size as user_row.
398The return value <code>hes[k]</code> is the second partial of
399\f$ w^{\rm T} F(x)\f$ with respect to the
400<code>row[k]</code> and <code>col[k]</code> component of \f$ x\f$.
401
402\param work
403This structure contains information that is computed by \c SparseHessianCompute.
404If the sparsity pattern, \c row vector, or \c col vectors
405are not the same between calls to \c SparseHessianCompute,
406\c work.clear() must be called to reinitialize \c work.
407
408\return
409Is the number of first order forward sweeps used to compute the
410requested Hessian values.
411(This is also equal to the number of second order reverse sweeps.)
412The total work, not counting the zero order
413forward sweep, or the time to combine computations, is proportional to this
414return value.
415*/
416template<class Base>
417template <class VectorBase, class VectorSet, class VectorSize>
418size_t ADFun<Base>::SparseHessianCompute(
419        const VectorBase&           x           ,
420        const VectorBase&           w           ,
421              VectorSet&            sparsity    ,
422        const VectorSize&           user_row    ,
423        const VectorSize&           user_col    ,
424              VectorBase&           hes         ,
425              sparse_hessian_work&  work        )
426{
427        using   CppAD::vectorBool;
428        size_t i, k, ell;
429
430        CppAD::vector<size_t>& row(work.row);
431        CppAD::vector<size_t>& col(work.col);
432        CppAD::vector<size_t>& color(work.color);
433        CppAD::vector<size_t>& order(work.order);
434
435        size_t n = Domain();
436
437        // some values
438        const Base zero(0);
439        const Base one(1);
440
441        // check VectorBase is Simple Vector class with Base type elements
442        CheckSimpleVector<Base, VectorBase>();
443
444        // number of components of Hessian that are required
445        size_t K = hes.size();
446        CPPAD_ASSERT_UNKNOWN( size_t( user_row.size() ) == K );
447        CPPAD_ASSERT_UNKNOWN( size_t( user_col.size() ) == K );
448
449        CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n );
450        CPPAD_ASSERT_UNKNOWN( color.size() == 0 || color.size() == n );
451        CPPAD_ASSERT_UNKNOWN( row.size() == 0   || row.size() == K );
452        CPPAD_ASSERT_UNKNOWN( col.size() == 0   || col.size() == K );
453
454
455        // Point at which we are evaluating the Hessian
456        Forward(0, x);
457
458        // check for case where nothing (except Forward above) to do
459        if( K == 0 )
460                return 0;
461
462        // Rows of the Hessian (i below) correspond to the forward mode index
463        // and columns (j below) correspond to the reverse mode index.
464        if( color.size() == 0 )
465        {
466                CPPAD_ASSERT_UNKNOWN( sparsity.n_set() ==  n );
467                CPPAD_ASSERT_UNKNOWN( sparsity.end() ==  n );
468
469                // copy user rwo and col to work space
470                row.resize(K);
471                col.resize(K);
472                for(k = 0; k < K; k++)
473                {       row[k] = user_row[k];
474                        col[k] = user_col[k];
475                }
476
477                // execute coloring algorithm
478                color.resize(n);
479                if( work.color_method == "cppad.general" )
480                        local::color_general_cppad(sparsity, row, col, color);
481                else if( work.color_method == "cppad.symmetric" )
482                        local::color_symmetric_cppad(sparsity, row, col, color);
483                else if( work.color_method == "colpack.general" )
484                {
485# if CPPAD_HAS_COLPACK
486                        local::color_general_colpack(sparsity, row, col, color);
487# else
488                        CPPAD_ASSERT_KNOWN(
489                                false,
490                                "SparseHessian: work.color_method = colpack.general "
491                                "and colpack_prefix missing from cmake command line."
492                        );
493# endif
494                }
495                else if(
496                        work.color_method == "colpack.symmetric" ||
497                        work.color_method == "colpack.star"
498                )
499                {
500# if CPPAD_HAS_COLPACK
501                        local::color_symmetric_colpack(sparsity, row, col, color);
502# else
503                        CPPAD_ASSERT_KNOWN(
504                                false,
505                                "SparseHessian: work.color_method is "
506                                "colpack.symmetric or colpack.star\n"
507                                "and colpack_prefix missing from cmake command line."
508                        );
509# endif
510                }
511                else
512                {       CPPAD_ASSERT_KNOWN(
513                                false,
514                                "SparseHessian: work.color_method is not valid."
515                        );
516                }
517
518                // put sorting indices in color order
519                VectorSize key(K);
520                order.resize(K);
521                for(k = 0; k < K; k++)
522                        key[k] = color[ row[k] ];
523                index_sort(key, order);
524
525        }
526        size_t n_color = 1;
527        for(ell = 0; ell < n; ell++) if( color[ell] < n )
528                n_color = std::max(n_color, color[ell] + 1);
529
530        // direction vector for calls to forward (rows of the Hessian)
531        VectorBase u(n);
532
533        // location for return values from reverse (columns of the Hessian)
534        VectorBase ddw(2 * n);
535
536        // initialize the return value
537        for(k = 0; k < K; k++)
538                hes[k] = zero;
539
540        // loop over colors
541# ifndef NDEBUG
542        const std::string& coloring = work.color_method;
543# endif
544        k = 0;
545        for(ell = 0; ell < n_color; ell++)
546        if( k == K )
547        {       // kludge because colpack returns colors that are not used
548                // (it does not know about the subset corresponding to row, col)
549                CPPAD_ASSERT_UNKNOWN(
550                        coloring == "colpack.general" ||
551                        coloring == "colpack.symmetic" ||
552                        coloring == "colpack.star"
553                );
554        }
555        else if( color[ row[ order[k] ] ] != ell )
556        {       // kludge because colpack returns colors that are not used
557                // (it does not know about the subset corresponding to row, col)
558                CPPAD_ASSERT_UNKNOWN(
559                        coloring == "colpack.general" ||
560                        coloring == "colpack.symmetic" ||
561                        coloring == "colpack.star"
562                );
563        }
564        else
565        {       CPPAD_ASSERT_UNKNOWN( color[ row[ order[k] ] ] == ell );
566
567                // combine all rows with this color
568                for(i = 0; i < n; i++)
569                {       u[i] = zero;
570                        if( color[i] == ell )
571                                u[i] = one;
572                }
573                // call forward mode for all these rows at once
574                Forward(1, u);
575
576                // evaluate derivative of w^T * F'(x) * u
577                ddw = Reverse(2, w);
578
579                // set the corresponding components of the result
580                while( k < K && color[ row[ order[k] ] ] == ell )
581                {       hes[ order[k] ] = ddw[ col[ order[k] ] * 2 + 1 ];
582                        k++;
583                }
584        }
585        return n_color;
586}
587// ===========================================================================
588// Public Member Functions
589// ===========================================================================
590/*!
591Compute user specified subset of a sparse Hessian.
592
593The C++ source code corresponding to this operation is
594\verbatim
595        SparceHessian(x, w, p, row, col, hes, work)
596\endverbatim
597
598\tparam Base
599is the base type for the recording that is stored in this ADFun<Base object.
600
601\tparam VectorBase
602is a simple vector class with elements of type \a Base.
603
604\tparam VectorSet
605is a simple vector class with elements of type
606\c bool or \c std::set<size_t>.
607
608\tparam VectorSize
609is a simple vector class with elements of type \c size_t.
610
611\param x [in]
612is a vector specifing the point at which to compute the Hessian.
613
614\param w [in]
615is the weighting vector that defines a scalar valued function by
616a weighted sum of the components of the vector valued function
617$latex F(x)$$.
618
619\param p [in]
620is the sparsity pattern for the Hessian that we are calculating.
621
622\param row [in]
623is the vector of row indices for the returned Hessian values.
624
625\param col [in]
626is the vector of columns indices for the returned Hessian values.
627It must have the same size are r.
628
629\param hes [out]
630is the vector of Hessian values.
631It must have the same size are r.
632The return value <code>hes[k]</code> is the second partial of
633\f$ w^{\rm T} F(x)\f$ with respect to the
634<code>row[k]</code> and <code>col[k]</code> component of \f$ x\f$.
635
636\param work
637This structure contains information that is computed by \c SparseHessianCompute.
638If the sparsity pattern, \c row vector, or \c col vectors
639are not the same between calls to \c SparseHessian,
640\c work.clear() must be called to reinitialize \c work.
641
642\return
643Is the number of first order forward sweeps used to compute the
644requested Hessian values.
645(This is also equal to the number of second order reverse sweeps.)
646The total work, not counting the zero order
647forward sweep, or the time to combine computations, is proportional to this
648return value.
649*/
650template<class Base>
651template <class VectorBase, class VectorSet, class VectorSize>
652size_t ADFun<Base>::SparseHessian(
653        const VectorBase&     x    ,
654        const VectorBase&     w    ,
655        const VectorSet&      p    ,
656        const VectorSize&     row  ,
657        const VectorSize&     col  ,
658        VectorBase&           hes  ,
659        sparse_hessian_work&  work )
660{
661        size_t n    = Domain();
662        size_t K    = hes.size();
663# ifndef NDEBUG
664        size_t k;
665        CPPAD_ASSERT_KNOWN(
666                size_t(x.size()) == n ,
667                "SparseHessian: size of x not equal domain dimension for f."
668        );
669        CPPAD_ASSERT_KNOWN(
670                size_t(row.size()) == K && size_t(col.size()) == K ,
671                "SparseHessian: either r or c does not have the same size as ehs."
672        );
673        CPPAD_ASSERT_KNOWN(
674                work.color.size() == 0 || work.color.size() == n,
675                "SparseHessian: invalid value in work."
676        );
677        for(k = 0; k < K; k++)
678        {       CPPAD_ASSERT_KNOWN(
679                        row[k] < n,
680                        "SparseHessian: invalid value in r."
681                );
682                CPPAD_ASSERT_KNOWN(
683                        col[k] < n,
684                        "SparseHessian: invalid value in c."
685                );
686        }
687        if( work.color.size() != 0 )
688                for(size_t j = 0; j < n; j++) CPPAD_ASSERT_KNOWN(
689                        work.color[j] <= n,
690                        "SparseHessian: invalid value in work."
691        );
692# endif
693        // check for case where there is nothing to compute
694        size_t n_sweep = 0;
695        if( K == 0 )
696                return n_sweep;
697
698        typedef typename VectorSet::value_type Set_type;
699        typedef typename local::internal_sparsity<Set_type>::pattern_type Pattern_type;
700        Pattern_type s;
701        if( work.color.size() == 0 )
702        {       bool transpose = false;
703                const char* error_msg = "SparseHessian: sparsity pattern"
704                " does not have proper row or column dimension";
705                sparsity_user2internal(s, p, n, n, transpose, error_msg);
706        }
707        n_sweep = SparseHessianCompute(x, w, s, row, col, hes, work);
708        return n_sweep;
709}
710/*!
711Compute a sparse Hessian.
712
713The C++ source code coresponding to this operation is
714\verbatim
715        hes = SparseHessian(x, w, p)
716\endverbatim
717
718
719\tparam Base
720is the base type for the recording that is stored in this
721ADFun<Base object.
722
723\tparam VectorBase
724is a simple vector class with elements of the \a Base.
725
726\tparam VectorSet
727is a simple vector class with elements of type
728\c bool or \c std::set<size_t>.
729
730\param x [in]
731is a vector specifing the point at which to compute the Hessian.
732
733\param w [in]
734The Hessian is computed for a weighted sum of the components
735of the function corresponding to this ADFun<Base> object.
736The argument \a w specifies the weights for each component.
737It must have size equal to the range dimension for this ADFun<Base> object.
738
739\param p [in]
740is a sparsity pattern for the Hessian.
741
742\return
743Will be a vector of size \c n * n containing the Hessian of
744at the point specified by \a x
745(where \c n is the domain dimension for this ADFun<Base> object).
746*/
747template <class Base>
748template <class VectorBase, class VectorSet>
749VectorBase ADFun<Base>::SparseHessian(
750        const VectorBase& x, const VectorBase& w, const VectorSet& p
751)
752{       size_t i, j, k;
753
754        size_t n = Domain();
755        VectorBase hes(n * n);
756
757        CPPAD_ASSERT_KNOWN(
758                size_t(x.size()) == n,
759                "SparseHessian: size of x not equal domain size for f."
760        );
761
762        typedef typename VectorSet::value_type Set_type;
763        typedef typename local::internal_sparsity<Set_type>::pattern_type Pattern_type;
764
765        // initialize the return value as zero
766        Base zero(0);
767        for(i = 0; i < n; i++)
768                for(j = 0; j < n; j++)
769                        hes[i * n + j] = zero;
770
771        // arguments to SparseHessianCompute
772        Pattern_type          s;
773        CppAD::vector<size_t> row;
774        CppAD::vector<size_t> col;
775        sparse_hessian_work   work;
776        bool transpose = false;
777        const char* error_msg = "SparseHessian: sparsity pattern"
778        " does not have proper row or column dimension";
779        sparsity_user2internal(s, p, n, n, transpose, error_msg);
780        k = 0;
781        for(i = 0; i < n; i++)
782        {       typename Pattern_type::const_iterator itr(s, i);
783                j = *itr;
784                while( j != s.end() )
785                {       row.push_back(i);
786                        col.push_back(j);
787                        k++;
788                        j = *(++itr);
789                }
790        }
791        size_t K = k;
792        VectorBase H(K);
793
794        // now we have folded this into the following case
795        SparseHessianCompute(x, w, s, row, col, H, work);
796
797        // now set the non-zero return values
798        for(k = 0; k < K; k++)
799                hes[ row[k] * n + col[k] ] = H[k];
800
801        return hes;
802}
803/*!
804Compute a sparse Hessian
805
806The C++ source code coresponding to this operation is
807\verbatim
808        hes = SparseHessian(x, w)
809\endverbatim
810
811
812\tparam Base
813is the base type for the recording that is stored in this
814ADFun<Base object.
815
816\tparam VectorBase
817is a simple vector class with elements of the \a Base.
818
819\param x [in]
820is a vector specifing the point at which to compute the Hessian.
821
822\param w [in]
823The Hessian is computed for a weighted sum of the components
824of the function corresponding to this ADFun<Base> object.
825The argument \a w specifies the weights for each component.
826It must have size equal to the range dimension for this ADFun<Base> object.
827
828\return
829Will be a vector of size \c n * n containing the Hessian of
830at the point specified by \a x
831(where \c n is the domain dimension for this ADFun<Base> object).
832*/
833template <class Base>
834template <class VectorBase>
835VectorBase ADFun<Base>::SparseHessian(const VectorBase &x, const VectorBase &w)
836{       size_t i, j, k;
837        typedef CppAD::vectorBool VectorBool;
838
839        size_t m = Range();
840        size_t n = Domain();
841
842        // determine the sparsity pattern p for Hessian of w^T F
843        VectorBool r(n * n);
844        for(j = 0; j < n; j++)
845        {       for(k = 0; k < n; k++)
846                        r[j * n + k] = false;
847                r[j * n + j] = true;
848        }
849        ForSparseJac(n, r);
850        //
851        VectorBool s(m);
852        for(i = 0; i < m; i++)
853                s[i] = w[i] != 0;
854        VectorBool p = RevSparseHes(n, s);
855
856        // compute sparse Hessian
857        return SparseHessian(x, w, p);
858}
859
860} // END_CPPAD_NAMESPACE
861# endif
Note: See TracBrowser for help on using the repository browser.