CppAD: A C++ Algorithmic Differentiation Package  20171217
sparse_hes_fun.hpp
Go to the documentation of this file.
1 # ifndef CPPAD_SPEED_SPARSE_HES_FUN_HPP
2 # define CPPAD_SPEED_SPARSE_HES_FUN_HPP
3
4 /* --------------------------------------------------------------------------
5 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
6
7 CppAD is distributed under multiple licenses. This distribution is under
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.
12 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
13 -------------------------------------------------------------------------- */
14 /*
15 $begin sparse_hes_fun$$16 spell 17 hes 18 cppad 19 hpp 20 fp 21 CppAD 22 namespace 23 const 24 bool 25 exp 26 arg 27$$ 28 29$section Evaluate a Function That Has a Sparse Hessian$$30 mindex sparse_hes_fun$$
31
32
33 $head Syntax$$34 codei%# include <cppad/speed/sparse_hes_fun.hpp> 35 %$$ 36$codei%sparse_hes_fun(%n%, %x%, %row%, %col%, %p%, %fp%)%$$37 38 head Purpose$$
39 This routine evaluates
40 $latex f(x)$$, latex f^{(1)} (x)$$, or$latex f^{(2)} (x)$$41 where the Hessian latex f^{(2)} (x)$$ is sparse.
42 The function $latex f : \B{R}^n \rightarrow \B{R}$$only depends on the 43 size and contents of the index vectors icode row$$ and$icode col$$. 44 The non-zero entries in the Hessian of this function have 45 one of the following forms: 46 latex $47 \DD{f}{x[row[k]]}{x[row[k]]} 48 \; , \; 49 \DD{f}{x[row[k]]}{x[col[k]]} 50 \; , \; 51 \DD{f}{x[col[k]]}{x[row[k]]} 52 \; , \; 53 \DD{f}{x[col[k]]}{x[col[k]]} 54$$$
55 for some $latex k $$between zero and latex K-1$$. 56 All the other terms of the Hessian are zero. 57 58$head Inclusion$$59 The template function code sparse_hes_fun$$
60 is defined in the $code CppAD$$namespace by including 61 the file code cppad/speed/sparse_hes_fun.hpp$$ 62 (relative to the CppAD distribution directory). 63 64$head Float$$65 The type icode Float$$ must be a $cref NumericType$$. 66 In addition, if icode y$$ and$icode z$$are icode Float$$ objects,
67 $codei% 68 %y% = exp(%z%) 69 %$$70 must set the icode y$$ equal the exponential of$icode z$$, i.e., 71 the derivative of icode y$$ with respect to $icode z$$is equal to icode y$$. 72 73$head FloatVector$$74 The type icode FloatVector$$ is any
75 $cref SimpleVector$$, or it can be a raw pointer, 76 with elements of type icode Float$$. 77 78$head n$$79 The argument icode n$$ has prototype
80 $codei% 81 size_t %n% 82 %$$83 It specifies the dimension for the domain space for latex f(x)$$. 84 85$head x$$86 The argument icode x$$ has prototype
87 $codei% 88 const %FloatVector%& %x% 89 %$$90 It contains the argument value for which the function, 91 or its derivative, is being evaluated. 92 We use latex n$$ to denote the size of the vector$icode x$$. 93 94 head row$$
95 The argument $icode row$$has prototype 96 codei% 97 const CppAD::vector<size_t>& %row% 98 %$$ 99 It specifies one of the first 100 index of$latex x$$for each non-zero Hessian term 101 (see cref/purpose/sparse_hes_fun/Purpose/$$ above).
102 All the elements of $icode row$$must be between zero and icode%n%-1%$$. 103 The value$latex K$$is defined by icode%K% = %row%.size()%$$.
104
105 $head col$$106 The argument icode col$$ has prototype 107$codei%
108  const CppAD::vector<size_t>& %col%
109 %$$110 and its size must be latex K$$; i.e., the same as for $icode col$$. 111 It specifies the second 112 index of latex x$$ for the non-zero Hessian terms. 113 All the elements of$icode col$$must be between zero and icode%n%-1%$$.
114 There are no duplicated entries requested, to be specific,
115 if $icode%k1% != %k2%$$then 116 codei% 117 ( %row%[%k1%] , %col%[%k1%] ) != ( %row%[%k2%] , %col%[%k2%] ) 118 %$$ 119 120$head p$$121 The argument icode p$$ has prototype
122 $codei% 123 size_t %p% 124 %$$125 It is either zero or two and 126 specifies the order of the derivative of latex f$$ 127 that is being evaluated, i.e.,$latex f^{(p)} (x)$$is evaluated. 128 129 head fp$$
130 The argument $icode fp$$has prototype 131 codei% 132 %FloatVector%& %fp% 133 %$$ 134 The input value of the elements of$icode fp$$does not matter. 135 136 subhead Function$$
137 If $icode p$$is zero, icode fp$$ has size one and 138$icode%fp%[0]%$$is the value of latex f(x)$$.
139
140 $subhead Hessian$$141 If icode p$$ is two,$icode fp$$has size icode K$$ and
142 for $latex k = 0 , \ldots , K-1$$, 143 latex $144 \DD{f}{ x[ \R{row}[k] ] }{ x[ \R{col}[k] ]} = fp [k] 145$$$ 146 147$children%
148  speed/example/sparse_hes_fun.cpp%
149  omh/sparse_hes_fun.omh
150 %$$151 152 head Example$$
153 The file
154 $cref sparse_hes_fun.cpp$$155 contains an example and test of code sparse_hes_fun.hpp$$. 156 It returns true if it succeeds and false otherwise. 157 158$head Source Code$$159 The file 160 cref sparse_hes_fun.hpp$$
161 contains the source code for this template function.
162
163 \$end
164 ------------------------------------------------------------------------------
165 */
166 // BEGIN C++
167 # include <cppad/core/cppad_assert.hpp>
169 # include <cppad/utility/vector.hpp>
170
171 // following needed by gcc under fedora 17 so that exp(double) is defined
172 # include <cppad/base_require.hpp>
173
174 namespace CppAD {
175  template <class Float, class FloatVector>
177  size_t n ,
178  const FloatVector& x ,
179  const CppAD::vector<size_t>& row ,
180  const CppAD::vector<size_t>& col ,
181  size_t p ,
182  FloatVector& fp )
183  {
184  // check numeric type specifications
185  CheckNumericType<Float>();
186
187  // check value of p
189  p == 0 || p == 2,
190  "sparse_hes_fun: p != 0 and p != 2"
191  );
192
193  size_t K = row.size();
194  size_t i, j, k;
195  if( p == 0 )
196  fp[0] = Float(0);
197  else
198  { for(k = 0; k < K; k++)
199  fp[k] = Float(0);
200  }
201
202  // determine which diagonal entries are present in row[k], col[k]
203  CppAD::vector<size_t> diagonal(n);
204  for(i = 0; i < n; i++)
205  diagonal[i] = K; // no diagonal entry for this row
206  for(k = 0; k < K; k++)
207  { if( row[k] == col[k] )
208  { CPPAD_ASSERT_UNKNOWN( diagonal[row[k]] == K );
209  // index of the diagonal entry
210  diagonal[ row[k] ] = k;
211  }
212  }
213
214  // determine which entries must be multiplied by a factor of two
215  CppAD::vector<Float> factor(K);
216  for(k = 0; k < K; k++)
217  { factor[k] = Float(1);
218  for(size_t k1 = 0; k1 < K; k1++)
219  { bool reflected = true;
220  reflected &= k != k1;
221  reflected &= row[k] != col[k];
222  reflected &= row[k] == col[k1];
223  reflected &= col[k] == row[k1];
224  if( reflected )
225  factor[k] = Float(2);
226  }
227  }
228
229  Float t;
230  for(k = 0; k < K; k++)
231  { i = row[k];
232  j = col[k];
233  t = exp( x[i] * x[j] );
234  switch(p)
235  {
236  case 0:
237  fp[0] += t;
238  break;
239
240  case 2:
241  if( i == j )
242  { // dt_dxi = 2.0 * xi * t
243  fp[k] += ( Float(2) + Float(4) * x[i] * x[i] ) * t;
244  }
245  else
246  { // dt_dxi = xj * t
247  fp[k] += factor[k] * ( Float(1) + x[i] * x[j] ) * t;
248  if( diagonal[i] != K )
249  { size_t ki = diagonal[i];
250  fp[ki] += x[j] * x[j] * t;
251  }
252  if( diagonal[j] != K )
253  { size_t kj = diagonal[j];
254  fp[kj] += x[i] * x[i] * t;
255  }
256  }
257  break;
258  }
259  }
260
261  }
262 }
263 // END C++
264 # endif
#define CPPAD_ASSERT_KNOWN(exp, msg)
Check that exp is true, if not print msg and terminate execution.
void sparse_hes_fun(size_t n, const FloatVector &x, const CppAD::vector< size_t > &row, const CppAD::vector< size_t > &col, size_t p, FloatVector &fp)
Define the CppAD error checking macros (all of which begin with CPPAD_ASSERT_)
AD< Base > exp(const AD< Base > &x)
size_t size(void) const
number of elements currently in this vector.
Definition: vector.hpp:387
#define CPPAD_ASSERT_UNKNOWN(exp)
Check that exp is true, if not terminate execution.
File used to define CppAD::vector and CppAD::vectorBool.