AdolcForward 4.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159
  1. // This file is part of Eigen, a lightweight C++ template library
  2. // for linear algebra.
  3. //
  4. // Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
  5. //
  6. // This Source Code Form is subject to the terms of the Mozilla
  7. // Public License v. 2.0. If a copy of the MPL was not distributed
  8. // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
  9. #ifndef EIGEN_ADLOC_FORWARD
  10. #define EIGEN_ADLOC_FORWARD
  11. //--------------------------------------------------------------------------------
  12. //
  13. // This file provides support for adolc's adouble type in forward mode.
  14. // ADOL-C is a C++ automatic differentiation library,
  15. // see https://projects.coin-or.org/ADOL-C for more information.
  16. //
  17. // Note that the maximal number of directions is controlled by
  18. // the preprocessor token NUMBER_DIRECTIONS. The default is 2.
  19. //
  20. //--------------------------------------------------------------------------------
  21. #define ADOLC_TAPELESS
  22. #ifndef NUMBER_DIRECTIONS
  23. # define NUMBER_DIRECTIONS 2
  24. #endif
  25. #include <adolc/adtl.h>
  26. // adolc defines some very stupid macros:
  27. #if defined(malloc)
  28. # undef malloc
  29. #endif
  30. #if defined(calloc)
  31. # undef calloc
  32. #endif
  33. #if defined(realloc)
  34. # undef realloc
  35. #endif
  36. #include "../../Eigen/Core"
  37. namespace Eigen {
  38. /**
  39. * \defgroup AdolcForward_Module Adolc forward module
  40. * This module provides support for adolc's adouble type in forward mode.
  41. * ADOL-C is a C++ automatic differentiation library,
  42. * see https://projects.coin-or.org/ADOL-C for more information.
  43. * It mainly consists in:
  44. * - a struct Eigen::NumTraits<adtl::adouble> specialization
  45. * - overloads of internal::* math function for adtl::adouble type.
  46. *
  47. * Note that the maximal number of directions is controlled by
  48. * the preprocessor token NUMBER_DIRECTIONS. The default is 2.
  49. *
  50. * \code
  51. * #include <unsupported/Eigen/AdolcSupport>
  52. * \endcode
  53. */
  54. //@{
  55. } // namespace Eigen
  56. // Eigen's require a few additional functions which must be defined in the same namespace
  57. // than the custom scalar type own namespace
  58. namespace adtl {
  59. inline const adouble& conj(const adouble& x) { return x; }
  60. inline const adouble& real(const adouble& x) { return x; }
  61. inline adouble imag(const adouble&) { return 0.; }
  62. inline adouble abs(const adouble& x) { return fabs(x); }
  63. inline adouble abs2(const adouble& x) { return x*x; }
  64. inline bool (isinf)(const adouble& x) { return (Eigen::numext::isinf)(x.getValue()); }
  65. inline bool (isnan)(const adouble& x) { return (Eigen::numext::isnan)(x.getValue()); }
  66. }
  67. namespace Eigen {
  68. template<> struct NumTraits<adtl::adouble>
  69. : NumTraits<double>
  70. {
  71. typedef adtl::adouble Real;
  72. typedef adtl::adouble NonInteger;
  73. typedef adtl::adouble Nested;
  74. enum {
  75. IsComplex = 0,
  76. IsInteger = 0,
  77. IsSigned = 1,
  78. RequireInitialization = 1,
  79. ReadCost = 1,
  80. AddCost = 1,
  81. MulCost = 1
  82. };
  83. };
  84. template<typename Functor> class AdolcForwardJacobian : public Functor
  85. {
  86. typedef adtl::adouble ActiveScalar;
  87. public:
  88. AdolcForwardJacobian() : Functor() {}
  89. AdolcForwardJacobian(const Functor& f) : Functor(f) {}
  90. // forward constructors
  91. template<typename T0>
  92. AdolcForwardJacobian(const T0& a0) : Functor(a0) {}
  93. template<typename T0, typename T1>
  94. AdolcForwardJacobian(const T0& a0, const T1& a1) : Functor(a0, a1) {}
  95. template<typename T0, typename T1, typename T2>
  96. AdolcForwardJacobian(const T0& a0, const T1& a1, const T1& a2) : Functor(a0, a1, a2) {}
  97. typedef typename Functor::InputType InputType;
  98. typedef typename Functor::ValueType ValueType;
  99. typedef typename Functor::JacobianType JacobianType;
  100. typedef Matrix<ActiveScalar, InputType::SizeAtCompileTime, 1> ActiveInput;
  101. typedef Matrix<ActiveScalar, ValueType::SizeAtCompileTime, 1> ActiveValue;
  102. void operator() (const InputType& x, ValueType* v, JacobianType* _jac) const
  103. {
  104. eigen_assert(v!=0);
  105. if (!_jac)
  106. {
  107. Functor::operator()(x, v);
  108. return;
  109. }
  110. JacobianType& jac = *_jac;
  111. ActiveInput ax = x.template cast<ActiveScalar>();
  112. ActiveValue av(jac.rows());
  113. for (int j=0; j<jac.cols(); j++)
  114. for (int i=0; i<jac.cols(); i++)
  115. ax[i].setADValue(j, i==j ? 1 : 0);
  116. Functor::operator()(ax, &av);
  117. for (int i=0; i<jac.rows(); i++)
  118. {
  119. (*v)[i] = av[i].getValue();
  120. for (int j=0; j<jac.cols(); j++)
  121. jac.coeffRef(i,j) = av[i].getADValue(j);
  122. }
  123. }
  124. protected:
  125. };
  126. //@}
  127. }
  128. #endif // EIGEN_ADLOC_FORWARD