Exotica
autodiff_chain_hessian_sparse.h
Go to the documentation of this file.
1 // Copyright (C) 2018
2 //
3 // This Source Code Form is subject to the terms of the Mozilla
4 // Public License v. 2.0. If a copy of the MPL was not distributed
5 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
6 //
7 // Modified from unsupported/Eigen/src/AutoDiff/AutoDiffJacobian.h
8 
9 #ifndef EIGEN_AUTODIFF_CHAIN_HESSIAN_SPARSE_H_
10 #define EIGEN_AUTODIFF_CHAIN_HESSIAN_SPARSE_H_
11 
15 
16 namespace Eigen
17 {
18 template <typename Functor>
20 {
21 public:
24 // forward constructors
25 #if EIGEN_HAS_VARIADIC_TEMPLATES
26  template <typename... T>
27  AutoDiffChainHessianSparse(const T &... Values) : Functor(Values...)
28  {
29  }
30 #else
31  template <typename T0>
32  AutoDiffChainHessianSparse(const T0 &a0) : Functor(a0)
33  {
34  }
35  template <typename T0, typename T1>
36  AutoDiffChainHessianSparse(const T0 &a0, const T1 &a1) : Functor(a0, a1)
37  {
38  }
39  template <typename T0, typename T1, typename T2>
40  AutoDiffChainHessianSparse(const T0 &a0, const T1 &a1, const T2 &a2) : Functor(a0, a1, a2)
41  {
42  }
43 #endif
44 
45  typedef typename Functor::InputType InputType;
46  typedef typename Functor::ValueType ValueType;
47  typedef typename ValueType::Scalar Scalar;
48 
49  enum
50  {
51  InputsAtCompileTime = InputType::RowsAtCompileTime,
52  ValuesAtCompileTime = ValueType::RowsAtCompileTime,
53  JacobianInputsAtCompileTime = Functor::JacobianColsAtCompileTime // JacobianInputsAtCompileTime no longer have to match InputsAtCompileTime
54  };
55 
56  typedef SparseMatrix<Scalar> JacobianType;
57 
58  typedef SparseMatrix<Scalar> InputJacobianType; // Jacobian.cols() matches InputJacobian.cols()
59  typedef Array<SparseMatrix<Scalar>, ValuesAtCompileTime, 1> HessianType;
60  typedef Array<SparseMatrix<Scalar>, InputsAtCompileTime, 1> InputHessianType;
61  typedef typename JacobianType::Index Index;
62 
63  typedef SparseVector<Scalar> InnerDerivativeType; // Derivative rows() matches InputJacobian.cols()
65  typedef SparseVector<InnerActiveScalar> OuterDerivativeType;
67 
68  typedef typename OuterDerivativeType::InnerIterator JacobianInnerIteratorType;
69  typedef typename InnerDerivativeType::InnerIterator HessianInnerIteratorType;
70 
71  typedef Matrix<OuterActiveScalar, InputsAtCompileTime, 1> ActiveInput;
72  typedef Matrix<OuterActiveScalar, ValuesAtCompileTime, 1> ActiveValue;
73 
74 #if EIGEN_HAS_VARIADIC_TEMPLATES
75  // Some compilers don't accept variadic parameters after a default parameter,
76  // i.e., we can't just write _jac=0 but we need to overload operator():
77  EIGEN_STRONG_INLINE
78  void operator()(const InputType &x, ValueType &v) const
79  {
80  this->operator()(x, v);
81  }
82 
83  template <typename... ParamsType>
84  void operator()(const InputType &x, ValueType &v, const ParamsType &... Params) const
85  {
86  this->operator()(x, v, Params...);
87  }
88 
89  template <typename... ParamsType>
90  void operator()(const InputType &x, ValueType &v, JacobianType &jac, const ParamsType &... Params) const
91  {
92  AutoDiffChainJacobianSparse<Functor> autoj(*static_cast<const Functor *>(this));
93  autoj(x, v, jac, Params...);
94  }
95 
96  template <typename... ParamsType>
97  void operator()(const InputType &x, ValueType &v, JacobianType &jac, const InputJacobianType &ijac,
98  const ParamsType &... Params) const
99  {
100  AutoDiffChainJacobianSparse<Functor> autoj(*static_cast<const Functor *>(this));
101  autoj(x, v, jac, ijac, Params...);
102  }
103 
104  template <typename... ParamsType>
105  void operator()(const InputType &x, ValueType &v, JacobianType &jac, HessianType &hess, const ParamsType &... Params) const
106  {
107  this->operator()(x, v, jac, hess, nullptr, nullptr, Params...);
108  }
109 
110  template <typename... ParamsType>
111  void operator()(const InputType &x, ValueType &v, JacobianType &jac, HessianType &hess, const InputJacobianType &ijac, const InputHessianType &ihess,
112  const ParamsType &... Params) const
113  {
114  this->operator()(x, v, jac, hess, &ijac, &ihess, Params...);
115  }
116 
117  // Optional parameter InputJacobian (_ijac)
118  template <typename... ParamsType>
119  void operator()(const InputType &x, ValueType &v, JacobianType &jac, HessianType &hess, const InputJacobianType *_ijac = 0, const InputHessianType *_ihess = 0,
120  const ParamsType &... Params) const
121 #else
122  EIGEN_STRONG_INLINE
123  void operator()(const InputType &x, ValueType &v) const
124  {
125  this->operator()(x, v);
126  }
127 
128  void operator()(const InputType &x, ValueType &v, JacobianType &jac) const
129  {
130  AutoDiffChainJacobianSparse<Functor> autoj(*static_cast<const Functor *>(this));
131  autoj(x, v, jac);
132  }
133 
134  void operator()(const InputType &x, ValueType &v, JacobianType &jac, const InputJacobianType &ijac) const
135  {
136  AutoDiffChainJacobianSparse<Functor> autoj(*static_cast<const Functor *>(this));
137  autoj(x, v, jac, ijac);
138  }
139 
140  void operator()(const InputType &x, ValueType &v, JacobianType &jac, HessianType &hess) const
141  {
142  this->operator()(x, v, jac, hess, nullptr, nullptr);
143  }
144 
145  void operator()(const InputType &x, ValueType &v, JacobianType &jac, HessianType &hess, const InputJacobianType &ijac, const InputHessianType &ihess) const
146  {
147  this->operator()(x, v, jac, hess, &ijac, &ihess);
148  }
149 
150  void operator()(const InputType &x, ValueType &v, JacobianType &jac = 0, HessianType &hess, const InputJacobianType *_ijac = 0, const InputHessianType *_ihess = 0) const
151 #endif
152  {
153  ActiveInput ax = x.template cast<OuterActiveScalar>();
154  ActiveValue av(jac.rows());
155 
156  // Provide either both input jacobian and hessian, or none
157  eigen_assert((_ijac && _ihess) || (!_ijac && !_ihess));
158 
159  if (!_ijac)
160  {
161  eigen_assert(x.rows() == jac.cols());
162 
163  for (Index j = 0; j < jac.rows(); ++j)
164  {
165  av[j].derivatives().resize(x.rows());
166  av[j].derivatives().reserve(x.rows());
167  for (Index k = 0; k < x.rows(); ++k)
168  av[j].derivatives().insert(k).derivatives().resize(x.rows());
169  }
170 
171  for (Index i = 0; i < x.rows(); ++i)
172  {
173  ax[i].derivatives().resize(x.rows());
174  ax[i].derivatives().insert(i) = 1.0;
175  ax[i].value().derivatives().resize(x.rows());
176  ax[i].value().derivatives().insert(i) = 1.0;
177  ax[i].derivatives().coeffRef(i).derivatives().resize(x.rows());
178  }
179  }
180  else
181  {
182  // If specified, copy derivatives from InputJacobian
183  const InputJacobianType &ijac = *_ijac;
184  const InputHessianType &ihess = *_ihess;
185 
186  eigen_assert(x.rows() == ihess.rows());
187  eigen_assert(ijac.cols() == ihess[0].rows() && ijac.cols() == ihess[0].cols());
188 
189  for (Index j = 0; j < jac.rows(); ++j)
190  {
191  av[j].derivatives().resize(ijac.cols());
192  av[j].derivatives().reserve(ijac.cols());
193  for (Index k = 0; k < ijac.cols(); ++k)
194  av[j].derivatives().insert(k).derivatives().resize(ijac.cols());
195  }
196 
197  for (Index i = 0; i < x.rows(); ++i)
198  {
199  ax[i].derivatives().resize(ijac.cols());
200  ax[i].derivatives() = ijac.row(i);
201  ax[i].value().derivatives().resize(ijac.cols());
202  ax[i].value().derivatives() = ijac.row(i);
203  for (Index k = 0; k < ijac.cols(); ++k)
204  {
205  ax[i].derivatives().coeffRef(k).derivatives() = ihess[i].row(k);
206  }
207  }
208  }
209 
210 #if EIGEN_HAS_VARIADIC_TEMPLATES
211  Functor::operator()(ax, av, Params...);
212 #else
213  Functor::operator()(ax, av);
214 #endif
215 
216  Index cols = _ijac ? _ijac->cols() : x.rows();
217  {
218  hess.resize(jac.rows());
219  for (Index i = 0; i < jac.rows(); ++i)
220  hess[i].resize(cols, cols);
221  }
222 
223  for (int i = 0; i < jac.rows(); ++i)
224  {
225  v[i] = av[i].value().value();
226  for (JacobianInnerIteratorType it(av[i].derivatives(), 0); it; ++it)
227  {
228  jac.insert(i, it.row()) = av[i].value().derivatives().coeffRef(it.row());
229  for (HessianInnerIteratorType ith(av[i].derivatives().coeffRef(it.row()).derivatives(), 0); ith; ++ith)
230  {
231  hess[i].insert(it.row(), ith.row()) = av[i].derivatives().coeffRef(it.row()).derivatives().coeffRef(ith.row());
232  }
233  }
234  }
235  }
236 };
237 
238 } // namespace Eigen
239 
240 #endif // EIGEN_AUTODIFF_CHAIN_HESSIAN_SPARSE_H_
autodiff_scalar.h
Eigen::AutoDiffChainHessianSparse::operator()
void operator()(const InputType &x, ValueType &v, JacobianType &jac, const InputJacobianType &ijac) const
Definition: autodiff_chain_hessian_sparse.h:134
Eigen
Definition: autodiff_chain_hessian.h:16
Eigen::AutoDiffChainHessianSparse::HessianType
Array< SparseMatrix< Scalar >, ValuesAtCompileTime, 1 > HessianType
Definition: autodiff_chain_hessian_sparse.h:59
Eigen::AutoDiffChainHessianSparse::Index
JacobianType::Index Index
Definition: autodiff_chain_hessian_sparse.h:61
Eigen::AutoDiffScalar
Definition: autodiff_scalar.h:48
Eigen::AutoDiffChainHessianSparse::InputsAtCompileTime
@ InputsAtCompileTime
Definition: autodiff_chain_hessian_sparse.h:51
Eigen::AutoDiffChainHessianSparse::ValueType
Functor::ValueType ValueType
Definition: autodiff_chain_hessian_sparse.h:46
exotica::Functor
FunctorBase< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::Dynamic > Functor
Definition: functor.h:49
Eigen::AutoDiffChainHessianSparse::AutoDiffChainHessianSparse
AutoDiffChainHessianSparse(const T0 &a0)
Definition: autodiff_chain_hessian_sparse.h:32
Eigen::AutoDiffChainHessianSparse::JacobianInnerIteratorType
OuterDerivativeType::InnerIterator JacobianInnerIteratorType
Definition: autodiff_chain_hessian_sparse.h:68
Eigen::AutoDiffChainHessianSparse::Scalar
ValueType::Scalar Scalar
Definition: autodiff_chain_hessian_sparse.h:47
Eigen::AutoDiffChainHessianSparse::operator()
EIGEN_STRONG_INLINE void operator()(const InputType &x, ValueType &v) const
Definition: autodiff_chain_hessian_sparse.h:123
Eigen::AutoDiffChainHessianSparse::InputType
Functor::InputType InputType
Definition: autodiff_chain_hessian_sparse.h:45
Eigen::AutoDiffChainHessianSparse::JacobianInputsAtCompileTime
@ JacobianInputsAtCompileTime
Definition: autodiff_chain_hessian_sparse.h:53
Eigen::AutoDiffChainHessianSparse::AutoDiffChainHessianSparse
AutoDiffChainHessianSparse(const T0 &a0, const T1 &a1)
Definition: autodiff_chain_hessian_sparse.h:36
Eigen::AutoDiffChainHessianSparse::JacobianType
SparseMatrix< Scalar > JacobianType
Definition: autodiff_chain_hessian_sparse.h:56
Eigen::AutoDiffChainHessianSparse::InnerActiveScalar
AutoDiffScalar< InnerDerivativeType > InnerActiveScalar
Definition: autodiff_chain_hessian_sparse.h:64
Eigen::AutoDiffChainHessianSparse::InnerDerivativeType
SparseVector< Scalar > InnerDerivativeType
Definition: autodiff_chain_hessian_sparse.h:63
Eigen::AutoDiffChainHessianSparse::ValuesAtCompileTime
@ ValuesAtCompileTime
Definition: autodiff_chain_hessian_sparse.h:52
Eigen::AutoDiffChainHessianSparse::OuterActiveScalar
AutoDiffScalar< OuterDerivativeType > OuterActiveScalar
Definition: autodiff_chain_hessian_sparse.h:66
Eigen::AutoDiffChainHessianSparse::ActiveValue
Matrix< OuterActiveScalar, ValuesAtCompileTime, 1 > ActiveValue
Definition: autodiff_chain_hessian_sparse.h:72
Eigen::AutoDiffChainHessianSparse::OuterDerivativeType
SparseVector< InnerActiveScalar > OuterDerivativeType
Definition: autodiff_chain_hessian_sparse.h:65
Eigen::AutoDiffChainHessianSparse::HessianInnerIteratorType
InnerDerivativeType::InnerIterator HessianInnerIteratorType
Definition: autodiff_chain_hessian_sparse.h:69
Eigen::AutoDiffChainHessianSparse::InputHessianType
Array< SparseMatrix< Scalar >, InputsAtCompileTime, 1 > InputHessianType
Definition: autodiff_chain_hessian_sparse.h:60
Eigen::AutoDiffChainJacobianSparse
Definition: autodiff_chain_jacobian_sparse.h:18
Eigen::AutoDiffChainHessianSparse::AutoDiffChainHessianSparse
AutoDiffChainHessianSparse(const T0 &a0, const T1 &a1, const T2 &a2)
Definition: autodiff_chain_hessian_sparse.h:40
Eigen::AutoDiffChainHessianSparse::operator()
void operator()(const InputType &x, ValueType &v, JacobianType &jac) const
Definition: autodiff_chain_hessian_sparse.h:128
Eigen::AutoDiffChainHessianSparse
Definition: autodiff_chain_hessian_sparse.h:19
Eigen::AutoDiffChainHessianSparse::operator()
void operator()(const InputType &x, ValueType &v, JacobianType &jac, HessianType &hess, const InputJacobianType &ijac, const InputHessianType &ihess) const
Definition: autodiff_chain_hessian_sparse.h:145
Eigen::AutoDiffChainHessianSparse::operator()
void operator()(const InputType &x, ValueType &v, JacobianType &jac, HessianType &hess) const
Definition: autodiff_chain_hessian_sparse.h:140
Eigen::AutoDiffChainHessianSparse::operator()
void operator()(const InputType &x, ValueType &v, JacobianType &jac=0, HessianType &hess, const InputJacobianType *_ijac=0, const InputHessianType *_ihess=0) const
Definition: autodiff_chain_hessian_sparse.h:150
Eigen::AutoDiffChainHessianSparse::InputJacobianType
SparseMatrix< Scalar > InputJacobianType
Definition: autodiff_chain_hessian_sparse.h:58
Eigen::AutoDiffChainHessianSparse::ActiveInput
Matrix< OuterActiveScalar, InputsAtCompileTime, 1 > ActiveInput
Definition: autodiff_chain_hessian_sparse.h:71
Eigen::AutoDiffChainHessianSparse::AutoDiffChainHessianSparse
AutoDiffChainHessianSparse(const Functor &f)
Definition: autodiff_chain_hessian_sparse.h:23
functor.h
Eigen::AutoDiffChainHessianSparse::AutoDiffChainHessianSparse
AutoDiffChainHessianSparse()
Definition: autodiff_chain_hessian_sparse.h:22
autodiff_chain_jacobian_sparse.h