Home | History | Annotate | Download | only in Eigen 
      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 (a] 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 
     10 #ifndef EIGEN_ADLOC_FORWARD
     11 #define EIGEN_ADLOC_FORWARD
     12 
     13 //--------------------------------------------------------------------------------
     14 //
     15 // This file provides support for adolc's adouble type in forward mode.
     16 // ADOL-C is a C++ automatic differentiation library,
     17 // see https://projects.coin-or.org/ADOL-C for more information.
     18 //
     19 // Note that the maximal number of directions is controlled by
     20 // the preprocessor token NUMBER_DIRECTIONS. The default is 2.
     21 //
     22 //--------------------------------------------------------------------------------
     23 
     24 #define ADOLC_TAPELESS
     25 #ifndef NUMBER_DIRECTIONS
     26 # define NUMBER_DIRECTIONS 2
     27 #endif
     28 #include <adolc/adouble.h>
     29 
     30 // adolc defines some very stupid macros:
     31 #if defined(malloc)
     32 # undef malloc
     33 #endif
     34 
     35 #if defined(calloc)
     36 # undef calloc
     37 #endif
     38 
     39 #if defined(realloc)
     40 # undef realloc
     41 #endif
     42 
     43 #include <Eigen/Core>
     44 
     45 namespace Eigen {
     46 
     47 /**
     48   * \defgroup AdolcForward_Module Adolc forward module
     49   * This module provides support for adolc's adouble type in forward mode.
     50   * ADOL-C is a C++ automatic differentiation library,
     51   * see https://projects.coin-or.org/ADOL-C for more information.
     52   * It mainly consists in:
     53   *  - a struct Eigen::NumTraits<adtl::adouble> specialization
     54   *  - overloads of internal::* math function for adtl::adouble type.
     55   *
     56   * Note that the maximal number of directions is controlled by
     57   * the preprocessor token NUMBER_DIRECTIONS. The default is 2.
     58   *
     59   * \code
     60   * #include <unsupported/Eigen/AdolcSupport>
     61   * \endcode
     62   */
     63   //@{
     64 
     65 } // namespace Eigen
     66 
     67 // Eigen's require a few additional functions which must be defined in the same namespace
     68 // than the custom scalar type own namespace
     69 namespace adtl {
     70 
     71 inline const adouble& conj(const adouble& x)  { return x; }
     72 inline const adouble& real(const adouble& x)  { return x; }
     73 inline adouble imag(const adouble&)    { return 0.; }
     74 inline adouble abs(const adouble&  x)  { return fabs(x); }
     75 inline adouble abs2(const adouble& x)  { return x*x; }
     76 
     77 }
     78 
     79 namespace Eigen {
     80 
     81 template<> struct NumTraits<adtl::adouble>
     82     : NumTraits<double>
     83 {
     84   typedef adtl::adouble Real;
     85   typedef adtl::adouble NonInteger;
     86   typedef adtl::adouble Nested;
     87   enum {
     88     IsComplex = 0,
     89     IsInteger = 0,
     90     IsSigned = 1,
     91     RequireInitialization = 1,
     92     ReadCost = 1,
     93     AddCost = 1,
     94     MulCost = 1
     95   };
     96 };
     97 
     98 template<typename Functor> class AdolcForwardJacobian : public Functor
     99 {
    100   typedef adtl::adouble ActiveScalar;
    101 public:
    102 
    103   AdolcForwardJacobian() : Functor() {}
    104   AdolcForwardJacobian(const Functor& f) : Functor(f) {}
    105 
    106   // forward constructors
    107   template<typename T0>
    108   AdolcForwardJacobian(const T0& a0) : Functor(a0) {}
    109   template<typename T0, typename T1>
    110   AdolcForwardJacobian(const T0& a0, const T1& a1) : Functor(a0, a1) {}
    111   template<typename T0, typename T1, typename T2>
    112   AdolcForwardJacobian(const T0& a0, const T1& a1, const T1& a2) : Functor(a0, a1, a2) {}
    113 
    114   typedef typename Functor::InputType InputType;
    115   typedef typename Functor::ValueType ValueType;
    116   typedef typename Functor::JacobianType JacobianType;
    117 
    118   typedef Matrix<ActiveScalar, InputType::SizeAtCompileTime, 1> ActiveInput;
    119   typedef Matrix<ActiveScalar, ValueType::SizeAtCompileTime, 1> ActiveValue;
    120 
    121   void operator() (const InputType& x, ValueType* v, JacobianType* _jac) const
    122   {
    123     eigen_assert(v!=0);
    124     if (!_jac)
    125     {
    126       Functor::operator()(x, v);
    127       return;
    128     }
    129 
    130     JacobianType& jac = *_jac;
    131 
    132     ActiveInput ax = x.template cast<ActiveScalar>();
    133     ActiveValue av(jac.rows());
    134 
    135     for (int j=0; j<jac.cols(); j++)
    136       for (int i=0; i<jac.cols(); i++)
    137         ax[i].setADValue(j, i==j ? 1 : 0);
    138 
    139     Functor::operator()(ax, &av);
    140 
    141     for (int i=0; i<jac.rows(); i++)
    142     {
    143       (*v)[i] = av[i].getValue();
    144       for (int j=0; j<jac.cols(); j++)
    145         jac.coeffRef(i,j) = av[i].getADValue(j);
    146     }
    147   }
    148 protected:
    149 
    150 };
    151 
    152 //@}
    153 
    154 }
    155 
    156 #endif // EIGEN_ADLOC_FORWARD
    157