cart-elc

Source code for CART-ELC
git clone git://git.laack.co/cart-elc.git
Log | Files | Refs | README | LICENSE

AdolcForward (4422B)

      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 
     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/adtl.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 inline bool (isinf)(const adouble& x) { return (Eigen::numext::isinf)(x.getValue()); }
     78 inline bool (isnan)(const adouble& x) { return (Eigen::numext::isnan)(x.getValue()); }
     79 
     80 }
     81 
     82 namespace Eigen {
     83 
     84 template<> struct NumTraits<adtl::adouble>
     85     : NumTraits<double>
     86 {
     87   typedef adtl::adouble Real;
     88   typedef adtl::adouble NonInteger;
     89   typedef adtl::adouble Nested;
     90   enum {
     91     IsComplex = 0,
     92     IsInteger = 0,
     93     IsSigned = 1,
     94     RequireInitialization = 1,
     95     ReadCost = 1,
     96     AddCost = 1,
     97     MulCost = 1
     98   };
     99 };
    100 
    101 template<typename Functor> class AdolcForwardJacobian : public Functor
    102 {
    103   typedef adtl::adouble ActiveScalar;
    104 public:
    105 
    106   AdolcForwardJacobian() : Functor() {}
    107   AdolcForwardJacobian(const Functor& f) : Functor(f) {}
    108 
    109   // forward constructors
    110   template<typename T0>
    111   AdolcForwardJacobian(const T0& a0) : Functor(a0) {}
    112   template<typename T0, typename T1>
    113   AdolcForwardJacobian(const T0& a0, const T1& a1) : Functor(a0, a1) {}
    114   template<typename T0, typename T1, typename T2>
    115   AdolcForwardJacobian(const T0& a0, const T1& a1, const T1& a2) : Functor(a0, a1, a2) {}
    116 
    117   typedef typename Functor::InputType InputType;
    118   typedef typename Functor::ValueType ValueType;
    119   typedef typename Functor::JacobianType JacobianType;
    120 
    121   typedef Matrix<ActiveScalar, InputType::SizeAtCompileTime, 1> ActiveInput;
    122   typedef Matrix<ActiveScalar, ValueType::SizeAtCompileTime, 1> ActiveValue;
    123 
    124   void operator() (const InputType& x, ValueType* v, JacobianType* _jac) const
    125   {
    126     eigen_assert(v!=0);
    127     if (!_jac)
    128     {
    129       Functor::operator()(x, v);
    130       return;
    131     }
    132 
    133     JacobianType& jac = *_jac;
    134 
    135     ActiveInput ax = x.template cast<ActiveScalar>();
    136     ActiveValue av(jac.rows());
    137 
    138     for (int j=0; j<jac.cols(); j++)
    139       for (int i=0; i<jac.cols(); i++)
    140         ax[i].setADValue(j, i==j ? 1 : 0);
    141 
    142     Functor::operator()(ax, &av);
    143 
    144     for (int i=0; i<jac.rows(); i++)
    145     {
    146       (*v)[i] = av[i].getValue();
    147       for (int j=0; j<jac.cols(); j++)
    148         jac.coeffRef(i,j) = av[i].getADValue(j);
    149     }
    150   }
    151 protected:
    152 
    153 };
    154 
    155 //@}
    156 
    157 }
    158 
    159 #endif // EIGEN_ADLOC_FORWARD